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I.  INTRODUCTION 


The  incomplete  beta  function  Ix(a,b)  is  defined  by 

Ix(a,b)  =G(a,b)rt^-'(l-t)*’-^dt,  a>0,  b>0,  0<x<l  (1) 

J  0 

B(a,b)=l/G(a,b)=f\^-‘(l-t)*»-‘dt=r(a)r(b)/r(a+b),  [10,  p.36],  (2) 

J  0 

where  the  gamma  function  r(u)  is  given  by 

fOO 

r(u)=  e“^  t“~^  dt,  u>0.  (3) 

J  0 


Thus  Ii(a,b)=l.  In  addition,  if  0<x<l  then  lx(0,b)  =  l  and  lx(a,0)  =  0. 


The  quantity  1  — Ix(a,b)  is  called  the  complement  of  Ix(a,b).  Using  u=l— t  in  (1)  and  (2) 

yields 

1-Ix(a,b)  =Iy(b,a),  y  =  l-x.  (4) 


The  function  Ix  occurs  in  many  branches  of  science,  including  atomic  physics,  fluid  dynamics, 
transmission  theory,  lattice  theory,  and  operations  research.  It  is  perhaps  best  known  for  its  extensive 
applications  in  statistics.  In  particular,  the  well-known  central  F-distribution  P(Folt^i,^2)  can  be 
obtained  from  Ix  by  the  following  substitutions  ; 


where 


a  =  (/i/2,  b  =  £/2/2  ,  t  =  i/,F/[«/2-t-«'iF]  ,  Fo= 


P{F,\u„.2)  K  +  r^,Fr^^'"^^)^^dF 


(5) 


(6) 


Q(f'ol  *^1'  *^2)  —  1  “  P(Pol  '^2)~Ji-x(*^2/2i  l'i/2). 


(7) 


The  incomplete  beta  function  is  also  directly  related  to  the  Student’s  t-distribution  A(to|i/)  and  the 
binomial  distribution  E(n,r,x),  where 

P(|t|<to)=A(tok)  =  ^  G(i,|)  (8) 

=  l-Ix(«//2,l/2),  ^  =  + 

E(n,r,x)=5];  ^P^x'Cl-x)"'*  =  Ix(r,n-r-(- 1).  (9) 

Derivations  of  these  well-known  results  are  given  in  [4]. 
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Procedures  for  computing  Ix  dale  back  to  Newton.  A  historical  survey  outlining  some  of  the 
analytical  and  approximation  methods  used  for  evaluating  Ix  is  given  by  Dutka  [7].  An  extensive 
literature  search  for  a  robust  algorithm  to  compute  Ix  did  not  reveal  a  publication  that  would  lead  to  a 
subroutine  acceptable  for  inclusion  in  a  high  quality  main  frame  mathematics  library  such  as  the 
NSWC  Library  of  Mathematics  Subroutines  (NSWCLIB),  [9].  The  best  procedure  found  was  in  a  report 
written  by  Amos  and  Daniel  [2].  For  a  special  case  however,  where  a  and  b  take  positive  integer  and 
half-integer  values,  an  algorithm  exists,  with  an  associated  Fortran  subroutine  ISUBX,  which  yields  Ix 
with  9-10  decimal-digit  accuracy,  [4,  5].  It  is  contained  in  NSWCLIB. 

In  this  report  an  algorithm  is  given  for  computing  lx  and  1—  Ix  -  A  transportable  Fortran 
subroutine  named  BRATIO  has  been  written  which  uses  the  algorithm.  BRATIO  is  designed  for  use  on 
computers  having  k-digit  single  precision  floating  arithmetics  where  6<k<14.  On  the  CDC  6000-7000 
series  computers,  BRATIO  yields  results  accurate  up  to  14  significant  digits  for  both  Ix  and  1— Ix- 
BRATIO  is  available  for  general  use  in  the  NSWC  mathematics  subroutine  library,  [9]. 

A  primary  region  of  difficulty  for  computing  Ix(a,b)  and  1— Ix(a,b)  has  been  when  a  and/or  b 
is  large  and  x«a/(a-(-b).  In  this  region  Ix  changes  rapidly  from  0  to  1.  A  continued  fraction,  with  new 
weighting  factors,  and  a  new  asymptotic  expansion  are  used  to  treat  this  region  satisfactorily. 

Section  II  contains  the  basic  equations  and  algorithms  used  for  BRATIO.  Section  III  describes 
in  a,  b,  x  space  the  regions  of  use  for  these  basic  relations.  An  associated  flowchart  for  BRATIO  is 
included.  Section  IV  describes  a  number  of  specialized  algorithms  required  in  order  to  use  the  basic 
relations  effectively.  Section  V  briefly  summarizes  the  accuracy  and  cfTiciency  of  BRATIO,  and  section 
VI  contains  a  few  examples  using  BRATIO.  Appendices  A,  B,  and  C  contain  proofs  for  some  results 
which  are  used  in  BRATIO.  A  Fortran  listing  of  BRATIO  and  its  required  subprograms  is  given  in 
Appendix  D. 
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II.  BASIC  RELATIONS 


Throughout  we  shall  use 

p=a/(a+b)  q  =  l-p=b/(a+b).  (10) 

Also  since  I^  and  1— Ix  arc  to  be  computed  to  the  greatest  possible  accuracy, 

y=l-x  (11) 

is  required  as  input  in  addition  to  a,  b,  and  x.  In  BRATIO,  relations  (12)  — (16)  are  used  to  compute 
either  lx  or  its  complement.  Their  domains  of  application  are  given  in  Section  III. 


BPSER 


c(a,b)«G(a,b)^(l+aE^- 


(j-b) 


J=> 


j!  (a+j ) 


(12) 


The  series  is  obtained  from  (1)  by  replacing  the  second  factor  in  the  integrand  with  its  binomial 
expansion.  Relation  (12)  is  used  only  when  x<0.7  and  b<l,  or  bx<0.7. 

Given  a  tolerance  f>0,  (12)  is  computed  by  the  function  BPSER(a,b,x,f ),  where 


Ix(a,b)«G(a,b)^(l+a;C^Wn  ) 

"'n  — Cn/(a  +  n),  Cn  =Ct).i(l  —  b/n)x,  Co=l 

N  =the  smallest  integer  such  that  a|wjyjl<e. 


(12.1) 

(12.2) 


BUP 


Ix(a,b)-Ix(a+  N,b)  =  x^'^  E  ~iy)r(t+j)^ 


(13) 


Given  a  tolerance  c,  (13)  is  computed  by  the  function  BUP(a,b,x,y,N,e).  BUP  is  used  only 
with  BPSER  or  the  next  relation  to  be  described,  BGRAT.  Relation  (13)  follows  from 


Ix(a+l,b)  =  Ix(a,b)  — x^y^G(a,b)/a,  which  can  be  obtained  by  substituting 
(a  +  b)t'‘{l-t)^-'  =  at''-‘(l-t)‘'-*-  ^[t"'  (1-t)'’] 


in 


Ix(a+  l,b)  — 


a+b 

aB(a,b) 


't^(l-t)*’-’dt. 


Equation  (13)  can  be  rewritten  in  the  form 

.a>  N-1 


(13.1) 


(13.2) 
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If  b<l  then  for  i>0,  dj_^^<dj  and  the  sequence  hj=djx'  is  monolonically  decreasing.  Thus,  the 

computation  of  the  sum  can  be  terminated  when  a  term  hm=dmX^  is  reached  that  satisfies 

m 

hm<  m=:N  —  1. 

i=o 


When  b>l  then  hj>hj^^  if  and  only  if  x<(a+l+i)/(a+b+i)  =  rj.  Also  we  note  that  rj<rj  when 
i<  j.  Thus  if  k  is  the  largest  integer  such  that  x>r|^_j,  then 

...  <h|^>h|^^j>  ...  >hf,_j. 

Therefore,  if  k  is  the  largest  integer  for  which  k<(b— l)x/y— a,  then  the  computation  of  the  sum  ^hj 

m 

can  be  terminated  when  a  term  h^  (m>k)  is  met  that  satisfies  h[p<  hp  of  m  =  N  — 1. 

i=:0 


BGRAT 


Ix(a,b)Rs  M  51!  Pn  Jn(b,u),  M  =  a>b 

n=o  r(a)T“ 


T=a+^-2^,  u=— Tlnx,  r=e~“  u“/r(b),  Jo(b,u)=  Q(b,u)/r 

n-1 


..-u  ..b , 


1  1 
Pn=(b— l)cn  +  n  l^^(nib  — n)  Cm  pn_m  ,  Cm  =  ^— 


.  Po  =  l 


Q(b,u)  = 


oo  -t  .b- 
c _ t 


V —  dt  (incomplete  gamma  function)  [1;  p.260] 
u  ^  (oJ 


Jn+i(b,u)  J„(b,u)+  (b^)'" 

L=smallest  integer  such  that  [  pj^Jj^l <  e  (  5Z  Pn  Jn +Wo/m)- 


(14) 

(14.1) 

(14.2) 

(14.3) 

(14.4) 

(14.5) 


Relation  (14)  is  used  only  when  a>15,  b<l,  and  x>0.7.  Given  e,  (14)  is  computed  by  the 
subroutine  BGRAT(a,b,x,y,w,(,lERR)  where  w  and  lERR  are  variables.  Given  an  initial  value  Wg  for 
w,  then  BGRAT  assigns  w  the  value  Wo  +  lx(a,b).  lERR  is  an  indicator  which  is  set  when  underflow 
forces  the  computation  of  (14)  to  end  prematurely. 


Equation  (14)  was  derived  by  Wise  [15].  Our  derivation  follows: 


-y/T 

If  t=e  ,  then 


‘(l-t)^-'dt  =  ^ 


b-i 


e  ^sinhb  i(X)dy 


(14.6) 


where  T  and  u  are  defined  at)ovc.  Substituting 
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sinh^-‘z- 


2j+b-i 


'  n=Q  '  j_T  J 

into  (H.6),  where  the  Cp  and  pj  satisfy  (14.2),  yields  (14)  where 

Jn(b,u)  =  ^^YpYj  Q(b+2n,u). 

Then  (14.4)  follows  from  (14.8)  and 

-u  b+2n  i  \ 

Q(b  +  2n  +  2,u)  “Q(b  +  2n,u)  + ^u  +  b+2n  +  l  j,  [1;  6.5.21] 
where  (14.9)  is  the  result  of  two  integrations  by  parts  of  its  left  hand  side. 

BFRAC 

)  B(a,b)  V  /?i+  /?,+  /Jm  -P“a+b 

Oi=l,  /li  =  ^^(A  +  l),  A=a-(a4-b)x  =  (a+b)(p-x) 

(a  +  n~l)(a  +  b4-n  —  1) 


(14.7) 


'n  +  i 


(a  +  2n-l)' 


■n(b— n)x‘‘,  n>l 


,  _  n(b-n)x  a  +  n 

n+)-"+a  +  2n-I  'a  +  2n+2L^^ 

H-n(l+y)]), 

n>0 

■n  +  i=  /^n  +  iAn  +  On  +  i  An-i 

Aq  =  0, 

A,  =a 

'n  +  i  —  ^n  +  i  f^n  +<*n  +  i  ^n-i 

cs 

o 

11 

Bi  =  /? 

rif]  U|  Ho  IJTfi 

Bn  0\'^  02  +  011 

m=  the  smallest  integer  such  that 

I  Am  /  Rrn  ~-^m-i  /^^m-il  ^  ^  l^m  /^inl- 


(14.8) 


(14.9) 


(15) 

(15.1) 

(15.2) 

(15.3) 

(15.4) 

(15.5) 


Equation  (15)  is  new  and  is  used  over  a  very  large  part  of  the  domain:  a>l,  b>40,  with  x<p. 
(liven  a  tolerance  (,  (15)  is  computed  by  the  function  BFRAC(a,b,x,y,A,().  The  weighting  factors  c„, 
given  below,  which  control  overflow  and  underflow,  appear  to  be  new. 

Relation  (15)  is  obtained  by  considering  the  cla-ssical  expansion 
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n  (b  — n) 


(a  +  2n-l)(a  +  2n) 


X,  n>0 


(15.6) 

(15.7) 


"a"-  131,  M  (,5.8) 

For  large  a  where  a;3>b,  we  note  that  djn^O  and  d2n  +  i '^  —  ( 1 +b/a)x.  Thus  if 
Rjj  is  the  iterate 

_  1  d,  dj  4n^2 

itn — i — T  i  T  ; — T 


*'"“!+  1+  1+  ■  ■  l+dn_,’ 

then  most  of  the  change  of  values  of  these  iterates  occurs  in  every  other  iterate.  Consequently,  the 

Af 


'associated”  continued  fraction  r-^  r— ^ 

b,+  boT 


a^  —  1 , 
b(  =  1  +  d,, 


is  considered  where 
^n  +  i  =  ~^2n-i  *^20 


^^n+i  =  >+d2n+d2n  +  i.  ">!•  [14;  P-20] 


The  iterates  are  the  iterates  for  this  expansion  as  already  pointed  out  by  Aroian  in  [3]  (  with  some 
corrections  given  in  [13]  ).  Hence 

Now  for  large  a  where  a;»b,  it  is  clear  that  d2n®s0  forces  anSsO.  However,  for  xssp 
we  also  find  that  d2n*5  0  forces  bn»sO.  This  can  cause  division  of  0  by  0  when  the  iterates  of  this 
continued  fraction  are  computed,  or  it  can  cause  division  to  overflow.  Thus,  this  continued  fraction  is 
also  not  considered  appropriate  for  computational  purposes. 

In  order  to  eliminate  the  problems  arising  from  d2n®^0.  we  rescale  the  coefficients  aj,  and  bj, 
with  weighting  factors  Cj,  such  that 


Oi— Cja,,  — c,|_|C|,af,  (n>2),  /7n— <^nl>n  (n>l)i 


where 


Cn=a-f2(n  — 1),  n>l. 

Then  the  iterates  of  ■  ■  ■  are  the  scaled  iterates  for  r— ^ 

0i+  02+  bl+  1^2  + 

obtain 


(15.9) 
and  we 


/  o,  o,  'I 

a  b 

1-  *  y  1 

r_L  -22_  ...  \ 

aB(a,b) 

vdj  +  02+  > 

B(a,b)  ' 

\0l  +02+  ' 

which  is  relation  (15).  The  expressions  (15. 1)  and  (15.5)  for  computing  (15)  follow  from  Theorem  2 

[1;  p.l9],  where  lim  -J—  ■■■  . 

^  ^  n-»oo  Bp  01+02  + 
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If  n<b,  then  and  is  a  positive  value  not  near  0.  To  insure  that  the  maximum 

number  of  iterations  satisfies  n<b,  (15)  is  applied  only  when  b>40.  In  this  case,  x  must  also  be  a 
sufficient  distance  from  p  when  a>  100. 


BASYM 


Let  7  be  an  arbitrary  positive  .scaling  factor  that  is  assigned  below.  Then 

Ix(a,b)«  e-z^  dn  Jn(z)  (/?7)".  x<P 
n=o 

b 


P  = 


a+b  ’ 


a  +  b’ 


0  =  . 


ab 


N(a  +  b)-^ 


U: 


a  b 

P  q 


B(a,b)N  ab 


2ir(a  +  b) 


(16)  I 

j 

I 


=  t,  t>0  (16.1) 

(yP(t)=a</)^l-,^  +  b<^^L^^,  0<t<l 

z  = 


a,!  —  7  "  A|i,  7>  0 

-'ll  =  [q  p  ■"  +  ( - 1 )"  p  q 

(Ao=l) 

(16.2) 

4'’=i 

(16.3) 

bj,^’  =  ran  +  ^  L 

i  =  i 

n=l,2,... 

Ill] 

{n>l) 

n-i 

do  =  I ,  dji  —  -  X!  <*i  ^„_i  + 1 
i=o 

(16.4) 

n_i 

=  cz- 

■oc 

e-v- v"  dv. 

V 

(16.5) 

l)-)  =:.J<|/a,  a<b 

/^7  =  Jp/b,  a>b. 

(16.6) 

N  is  defined  as  the  smallest  integer  such  that 

N 

^  II  mO 


i 


Relation  (16)  is  used  only  when  a  and  b  are  large  and  x«p.  It  appears  to  be  new.  Given 
A  =  (a  +  b)(p  — x)  and  a  tolerance  c,  (16)  is  computed  by  the  function  BASYM(a,b,x,y,A,c).  The 
expansion  is  obtained  by  writing  (1)  in  the  form 


a  b 

lx(a,b)=£J- 


’B(a,b) 


f’'  (l-t\^  dt 

Jo  I  ‘I  i  t(l- 


t)’ 


0<x<l. 


(16.7) 


From  (16.1),  v?(t)>0  for  t^p  and 

'^(x)  =- (a  In^  +  b  In^^  (16.8) 

Thus  if  u  =  J^(t)  for  t<p  and  u  =  — ,  ^(t)  for  t>p,  then 


Ix(a,b)  =  U/? 


du, 


(16.9) 


where  z=  Ji^(x)  if  x<p  and  z  =  —  ^^^(x)  if  x>p.  Hereafter,  we  shall  assume  that  x<p.  In  order  to  use 
(16.9),  u/(p  — t)  will  be  replaced  by  its  Maclaurin  series  in  u. 

From  (16.8)  and 


a  lnp=a  ln( 

P-t^ 

n=i 

,  P  / 

1  P  J 

b  In  =b  In  ( 

'l  +  Eni' 

,  oo  , 

)=-^Zh 

'  n=i 

one  obtains 

u2=^(t)  =  ;^(p-t)-5]  An(p-t)",  )p-t|<min{p,q}, 

n>o 

where  Ap  is  given  in  (16.2).  Thus  if  an  =  7“"An  for  7>0,  then  2(/37u)^=s^A,  where 

A  =  5I^ns”  s  =  7(P~t)'  Hence,  if  P(s)=.s>rA  then  T2/?7U  =  P(s).  If  v  =  P(s),  since 

n>o 

P  is  analytic  at  0  with  P(0)  =  0,  by  the  Lagrange-Biirmann  expansion  [8;  p.58]  the  inverse 

s  =  P“'(v)  of  P  is  given  by 

_ _ ,  ,,n 

s  —  2^  Cn  V  , 
n>i 

«'n=nrc8(P'‘‘) 

where  resfP”'')  is  the  residue  of  the  series  P”*'.  Also  for  any  r^tO 

kX) 
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where  is  given  in  (16.3),  [11].  Consequently, 


_1  (-n/2) 
^n  —  n  "n-i 


p  ^  n>i 


n>o 


where  dn  is  given  in  (16. d),  [11]. 

Now  if  h=a/b  for  a<b  and  h  =  b/a  for  a>b,  let  7  =  (l+h)/h,  (see  (16.6)).  Then  A  =  ^'ans’^ 

n>o 

for  |s|  <  1  and 


»n  pO  “  * )"  +  »>b. 


(16.10) 


It  can  be  shown  for  a  =  b  tliat  tlio  series  v/s  =  52  <lnv”  has  a  radius  of  convergence  no  larger  than  >l27r. 

n>o 

Indeed,  we  have 

=av*’ =  — a  ln(  1  — s*), 

and,  with  v  complex,  s  =  0  when  v  =  '(ix  exp(ix’/4). 

If  a<b  then  A  =  a„.s'’>0  for  0<.s<l  since  each  an>0.  Hence  vt=S'iA  is  defined  and  real  for 

n>o 

0<s<l.  Also,  from  (16.8)  we  note  that  v-*0C'  whens-»l,  and  dv/ds  =  (s/v)[(  1  — s)(  l+hs)]"' ^0  for 
0<s<l.  Consequently,  f(v)  =  v7s  =  >rA  can  be  regarded  as  a  function  of  v  for  all  v>0.  Since  it  is  also 
true  that  the  derivatives  f"''(v)  are  bounded  for  v>0  (see  Appendix  C),  let  =sup{  lf*”*(v)l :  v>0) 
for  n>  I.  Then  M  i<  1  (.see  Apiiendix  C),  and  for 

t''n{v)  =  f(v)-^  djv‘. 

i=o 


|i/'n(v)|<Mn  v"/n!  (v>0)  by  the  Taylor  formula  with  remainder.  Therefore,  from  ■f2/?u/(p  — t)=f(v) 
and  (16.9)  we  obtain 


I;,(a,b)  =  U^e- 


'£'d;(/f7)7li(2.)  +Kn 
i=o 


z>0 


*■’11  =  .^'  "  rii(>f2/f7u)dii, 

where  Jn(*)  's  given  by  (16.r)).  Also 

|Hnl<Mn/"!(^^7)'’Jn(2)-Mn/n!Jn(*)Nl+l>)]""'^'  (16.11) 
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so  that  (En|-*0  when  a-»oo  for  fixed  h.  Consequently,  (16)  is  asymptotic  for  a<b.  Also  we  note  that 
Jn(z)=2‘"''"'“V'r[(n  +  l)/2]  Q[(n  +  l)/2,z2)l, 
where  Q  is  defined  by  (14.3),  and  that  J|,(z)  can  be  computed  recursively  by 
Jo(z)  =  (^ri/4)e*^rfc(z),  = 

Jn(z)=2-^/'(,r2z)’'-'  +  (n-l)Jn_2(z),  n>2. 

When  a>b  the  situation  is  less  satisfactory,  since  it  has  not  yet  been  shown  that  (16)  is 
asymptotic.  Nevertheless,  the  utility  of  (16)  has  been  established  for  y<1.05q  when  b>100  by  extensive 
computer  testing. 

Finally,  we  observe  that  the  defintion  of  U  given  in  (16)  is  not  suitable  for  computational 
purposes.  U  can  be  accurately  evaluated  using 

A(a)  =ln  r(a)  — (a— |)ln  a+a— jln  (2)r)  (16.12) 

In  U=A(a  +  b)- A(a)-A(b). 


10 


111.  DOMAINS  FOR  SUBPROGRAMS 


In  this  section  we  specify  the  regions  of  application  for  the  five  subprograms  discussed  in  the 
previous  section.  A  flowchart  is  given  in  Fig.  1. 

In  order  to  establish  such  regions  the  following  conditions  must  be  met: 

(a)  The  applied  algorithms  must  be  efilcicnt  and  yield 
the  desired  accuracy  over  such  regions. 

(b)  It  is  necessary  that  J<.9,  where  J  has  the  value 
lx(a,b)  or  ly(b,a).  J  is  computed  by  the  proper 
choice,  acccording  to  (a),  of  one  or  more  of  the 
subprograms  based  on  (12)  — (16).  The  complement 
is  then  obtained  as  1— J. 

Even  though  analysis  was  carried  out  to  predict  the  efficiency  and  achievable  accuracy  of  the 
basic  relations  over  various  domains,  (a)  was  established  with  exhaustive  testing  by  Morris  using 
double  precision  versions  of  the  subprograms. 

Since  J  is  to  be  computed  to  the  greatest  possible  accuracy,  the  relative  tolerance  to  be  satisfied 
will  be  c  =  max{eo,10“'^}  where  Cq  **  smallest  number  for  which  ^he  floating  point 

arithmetic  being  used.  The  restriction  that  <>10”’^,  thereby  limiting  the  maximum  precision  to 
14  —  15  significant  digits,  is  made  since  many  of  the  supporting  subprograms  are  accurate  to  a 
maximum  of  14  digits  (see  section  IV). 

Arguments  are  presented  near  the  end  of  this  section  showing  that  condition  (b)  is  always 
satisfied.  For  ea.sy  reference  we  have: 

BPSER(a,b,x,<')  Relalion(12) 

BUP(a,b,x,y,N,f )  Relation(13) 

BGRAT(a,b,x,y,w.l5(,lERR)  Relation(  14) 

BFRAC(a,b,x,y,A,1.5f )  Relalion(  15) 

BASYM(a,b,x,y,A,l()(b ) '  Rclation(  16) 

There  are  two  main  domains  to  consider,  namely  min(a,b)<,l  and  min(a,b)>l.  It  should  be 
recalled  that  if  a  and  b.  and  x  and  y  are  interchanged  in  one  of  the  above  subprograms,  then  the 
subprogram,  except  for  BUP,  yields  J  =  Iy(b,a)=  1  — !x(a,b)  with  1— J=lx(a,b).  BUP  gives  on  the 
interchange  Iy(b,a)  —  Iy(b+  N,a). 


II 


min(a,b)<l 


If  x>l/2,  then  a  and  b,  and  x  and  y  are  interchanged.  For  x<l/2,  (17)  is  used  for  computing 
J=Ix(a,b)  and  (18)  — (20)  are  used  for  computing  J=ly(a,b).  In  (19)  and  (20),  Wq  is  the  initial  value 
of  w  and  J  is  the  final  value  of  w. 


BPSER(a,b,x,€)  max(a,b)<l,  a>min(0.2,b)  (17) 

max(a,b)<l,  a<min(0.2,b),  x^<0.9  (17.1) 

max(a,b)>l,  b<l  (17-2) 

max(a,b)>l,  b>l,  x<0.1,  (bx)^<0.7  (17.3) 

BPSER(b,a,y,£)  max(a,b)<l,  a<min(0.2,b),  x^>0.9,  x>0.3  (18) 

max(a,b)>l,  b>l,  x>0.3  (18.1) 

BGRAT(b,a,y,x,w,15£,IERR),  wo=0 

max(a,b)>l,  b>15,  0.1<x<0.3  (19) 

max(a,b)>l,  b>15,  x<0.1,  (bx)^>0.7  (191) 

BGRAT(b  +  N,a,y,x,w,15€,IERR),  Wo  =  BUP(b,a,y,x,N,£),  N  =  20 

max(a,b)>l,  b>l,  0.1<x<0.3,  b<15  (20) 

max(a,b)>l,  b>I,  x<0.1,  (bx)^>0.7,  b<15  (20.1) 

max(a,b)<l,  a<min(0.2,b),  x*>0.9,  x<0.3  (20.2) 


min(a,b)>l 

If  x>p  then  a  and  b,  and  x  and  y  are  interchanged.  For  x<p,  (21)  — (26)  are  used  for 
computing  J  =  Ix(a,b).  In  (22)  — (24)  N  is  the  largest  integer  less  than  b  and  b  =  b  — N.  Also,  in  (23)  and 
(24)  Wg  is  the  initial  value  of  w  and  J  is  the  final  value  of  w. 

BPSER(a,b,x,f)  b<40,  bx<.7  (21) 

BUP(b,a,y,x,N,f)  +  BPSER(a,b,x,f) 

b<40,  bx>.7,  x<.7  (22) 

BGRAT(a,b,x,y,w,15f  ,IERR),  Wo  =  BUP(b,a,y,x,N,() 

b<40,  x>.7,  a>15  (23) 
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BGRAT(a+M,b,x,y,w,15f,IERR),  M=20 
Wo  =  BUP(b,a,y,x,N,f)  +  BUP(a,b,x,y,M,() 


b<40. 

x>.7. 

a<15 

(24) 

BASYM(a,b,x,y,A,100c) 

b>40. 

100<a<b, 

x>.97p 

(25) 

b>40. 

100<b<a, 

y  <1.03q 

(25.1) 

BFRAC(a,b,x,y,A,15f) 

b>40, 

a<b. 

a<100 

(26) 

b>40, 

100<a<b, 

x<.97p 

(26.1) 

b>40. 

a>b. 

b<100 

(26.2) 

b>40. 

100<b<a, 

y>1.03q 

(26.3) 

These  statements  arc  summarized  in  the  flowchart  for  BRATIO  in  Fig.l.  Proofs  are  now  given 
which  verify  that  J<0.9  is  always  satisfied  (requirement  b  above).  The  arguments  use  the  facts  that 
Inr(t)  is  strictly  convex  for  t>0,  [1;  6.4.10],  and  that  Ix(a,b)  is  a  decreasing  function  of  a  and  an 
increasing  function  of  b.  The  latter  result  is  proven  in  Appendix  B.  Since  In  r(t)  is  strictly  convex  for 
t>0,  we  also  note  that  ti’(t)=^{ln  r(t)},  (1;  6.3.1],  is  an  increasing  function  of  t. 

For  (17):  a<l,  b<l,  x<l/2,  a>min(.2,b) 

If  a>0.2  then 

J=Ix(a,b)<Ix(a,l)  =  x'*<(l/2)-^  =  0.8706. 

If  a>b  then 

J  =  Ix(aib)<Ix(a,a)<Ij^^(a,a). 

Also  I^^^(a,a)  =  0.5  from  (4). 

For  (17.1):  a<l,  b<l,  x<l/2,  a<min(.2,b),  x®<0.9 
J=Ix(a,b)<Ix(a,l)  =  x*<0.9 
For  (17.2):  a>I,  b<l,  x<l/2 

J  =  Ix(a,b)<Ix(l,l)  =  x<l/2 

For  (17.3):  a<l,  b>l,  x<0.1,  (bx)*<0.7 
Integrating  (1)  by  parts  gives 

Ix(a,b)=C(a,b)|^(  1  -x)^^-'  I*  t«(l  -t)'’-^dt| 

<C:(a,b)|^  (1  -x)'’-  '  x'‘|^  (1  -t)’’-"dt|=  G(a,b)^. 
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F(a,b)=r(a+b)/[r(b)b^],  0<a<l,  b>l. 

Since  5/5a[lnF]  =  ^(a+b)  —  lnb  is  an  increasing  function  of  a  for  a>0,  InF  is  strictly  convex  and 
hence  F  is  strictly  convex  for  a>0.  Thus  F(a,b)<l  for  0<a<l  since  F(0,b)  =  F(l,b)  =  l.  Hence 
J=Ix(a,b)<(bxf /r(a+l)<0.7/r(1.46163-)<0.791. 

For  (18):  a<l,  b<l,  0.3<x<0.5,  a<min(.2,b),  x^>0.9 
From  (1),  with 

(l-t)*’“'>l,  0<t<x, 

Ix(a,b)>H(a)^x^>H(a)(0.9)/2>.45,  (a<b), 

where 

H(a)  =r(a+b  +  l)/[r(a+l)r(b+l)]  (0<a<l). 

The  last  inequality  on  Ix  follows  from  the  fact  that  H(a)  is  increasing  and  that  H(0)  =  1.  Hence 

J  =Iy(b,a)  =  l-Ix(a,b)<0.55. 

For  (18.1):  a<l,  b>l,  0.3<x<0.5 

J=Iy(b,a)<Iy(l,l)  =  l-x<0.7 


For  (19):  a<l,  b>15,  0.1<x<0.3 

J=Iy(b,a)<Iy(b,l)=y^  =  (l-x)'’<(0.9)‘®«0.2059 


For  (19.1):  a<l,  b>15,  x<0.1,  (bx)®>0.7 


Then 


J=Iy(b,a)<Iy(b,l)  =  y^=(l-x)‘’. 


From(l),  using  (1—t)^  '>(1— x)^  * 


Ix(a,b)  > 


F(a,b)  (1-x)^ 
r(a)  1-x 


and  recalling  F(a,b)  in  the  proof  of  (17.3), 
F(a,b)(bx)^J/r(a+l). 


Now  5F/9b=g(a,b)F  where  g(a,b)  =  V>(a+b)  — ^(b)— a/b.  Since  gaa<0  (1;  6.4.10],  g  is  strictly  concave 
for  0<a<l.  Thus,  since  g(0,b)=g(l,b)=0  [1;  6.3.5],  g(a,b)>0  for  0<a<l  and  5F/5b>0.  Hence,  F  is 
increasing  in  b  for  0<a<l  and 


F(a,b)>F(a,l)  =  r(a+l). 


Thus  Ix(a,b)>0.7J  for  0<a<l,  which  implies  that  J<l/l.7ss.59. 


For  (20):  a<l,  l<b<15,  0.1<x<0.3 

J=Iy(b,a)<Iy(b,l)=y^=(l-x)^<(0.9)’^<0.9 

For  (20.1):  a<l,  l<b<15,  x<0.1,  (bx)®>0.7 
Proof  is  the  same  as  that  given  for  (19.1). 

For  (20.2):  b<l,  x<0.3,  a<min(.2,  b),  x^>0.9 

Ix(a,b)>Ix(a,a)  =  K(a)  ^  t^'  ‘  (1  -t)^'"  dt, 

where 

K(a)  =r(2a+l)/[r(a+l)r(a+l)],  (0<a<l). 

From  (i),  using 

(l-t)'‘"'>l,  0<t<x, 

it  follows  that 

Ix(a,b)>K(a)x®/2>0.45K(a). 

Since  <9K(a)/3a=2K(a)[t//(2a  +  l)  — ti>(a4-l)]>0,  K(a)  is  an  increasing  function.  Hence  K(a)>l  and 

J  =Iy(b,a)<0.55. 

For  (21)-(26):  a>l,  b>l 
In  this  case  it  can  be  shown  that  if  x<p  then 

fl  — 1/e  a<b 

Ix(a,b) <  ^  (27) 

(1/2  a>b. 

A  proof  of  this  result  is  given  in  Appendix  A.  Hence 

J=  f>x(a.b)  x<p 

\ly(b,a)  x>p. 
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Notation  for  the  flowchart  in  Fig.l 


12  Refers  to  (12)  for  computing  lx(a,b)  (BPSER) 

12  Refers  to  (12)  for  computing  Iy(b,a)  (BPSER) 

13  Refers  to  (13)  for  computing  Ix(a,b)  (BUP) 

13  Refers  to  (13)  for  computing  Iy(b,a)  (BUP) 

14  Refers  to  (14)  for  computing  Ix(a,b)  (BGRAT) 

14  Refers  to  (14)  for  computing  Iy(b,a)  (BGRAT) 

15  Refers  to  (15)  for  computing  Ix(a,b)  (BFRAC) 

15  Refers  to  (15)  for  computing  Iy(b,a)  (BFRAC) 

16  Refers  to  (16)  for  computing  Ix(a,b)  (BASYM) 

16  Refers  to  (16)  for  computing  Iy(b,a)  (BASYM) 

p=a/(a4-b),  q  =  b/(a+b) 

I.C.^  Interchange  a  and  b;  x  and  y. 

[b]=largest  integer  <  b. 

A=a— (a+b)x,  a<b;  A  =  (a4-b)y  — b,  a>b. 

Numbers  above  some  of  the  flowchart  boxes  refer  to  the  labels  in  the  Fortran  listing  of 
BRATIO. 
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Flowchart  for  BRATIO  (Notation 


IV.  AUXILIARY  FUNCTIONS 


In  order  to  compute  Ix(a,b)  and  1  — Ix(a,b),  procedures  are  needed  for  evaluating  r(a),  In  r(a), 
the  error  function  erf  x,  exp(x^)erfc  x,  the  incomplete  gamma  function  Q(a,x),  (14.3),  for  a<l,  and  the 
functions 


e^-1 


ln(l+x) 

In  r(l  4-x) 

l/r(l+a)-l 

c-’‘x®/r(a) 

^(x)  =  x—  1  — Inx 


(29) 


(|x|<.375) 
(-0.2<a<1.25) 
(-0.5<a<1.5) 
(a>0,  x>0) 
(x>0). 


These  functions  arc  discussed  in  [6].  Also,  procedures  are  needed  for  computing 


A(a)  =  ln  r(a)  — (a  — .5)ln  a+a— .5  ln(27r) 
ALGDlV(a,b)  =  ln[r(b)/r(a+b)l, 
BCORR(a,b)  =  A(a)  +  A(b)  -  A(a+ b) 
BETALN(a,b)  =  ln  B(a,b) 
BRCOMP(a,b,x,y )  =x^y‘^  /  B(a,b) 


(a>8) 

(a>0,b>8) 

(a,b>8) 

(a,b>0) 

(a,b>0,  0<x<l,  y  =  l— x). 


(30) 


Rational  minimax  approximations  are  used  for  the  functions  given  in  (29).  Experience  indicates  that 
such  approximations  normally  generate  less  error  and  can  be  considerably  more  efficient  than  the 
standard  expansions.  However,  minimax  approximations  have  the  disadvantage  of  being  limited  to  a 
fixed  maximum  precision.  The  mimimax  approximations  used  are  designed  to  achieve  a  maximum 
precision  of  14  significant  digits. 


If  A(a)  is  needed  only  for  a >20,  then  the  sum 
l/(12a)-l/(360a^)+ 

in  the  asymptotic  expansion  of  In  r(a)  [1;  6.1.41]  may  be  used.  If  a>15  then  the  minimax 
approximation 

A(a)=  X:  c„/a^"-^‘ 
n  =  o 


Co=  .83333  33333  33.333E-01 
c,=  -. 27777  77777  70481E-02 
C2=  .79365  06631  83693E-03 
C3=-..59515  63.364  28591E-03 
c,,=  .82075  63703  5.3826E-03 
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can  be  applied,  and  if  a>8  the  minimax  approximation 

A(a)=  E  dn/a'"^' 
n  =  o 


do=  .83333  33333  33333E-01  (32) 

di= -.27777  77777  60991E-02 

d2=  .79365  06668  25390E-03 

d3= -.59520  29313  51870E-03 

d^=  .83730  80340  312I5E-03 

d5=-. 16532  29627  80713E-02 


can  be  used.  These  approximations  were  obtained  by  Morris  [9].  On  the  CDC  6000-7000  series 
computers,  they  are  accurate  to  within  1  unit  of  the  14‘*  significant  digit. 

Expansions  for  ALGDIV(a,b)  and  BCORR(a,b)  use  (32).  From  the  definition  of  A 

ALGDIV(a,b)=w  — (a  +  b  — .5)ln(l +a/b)— a(ln  b— 1)  (33) 

w=A(b)-A(a  +  b). 

Let 

p=a/(a-f-b),  q  =  b/(a-f-b),  Sm  =  1 +q+ •  ■  ■ +q^  *  (m>l). 

Then 

l-q'"=(l-,)S„=pS„„  ^  =  ,3,) 

Thus,  from  (32)  we  obtain 


5 

E  ‘^n 


n=o 


^2n+i 
b2n  > 


(35) 


which  completes  the  algorithm  for  ALGDIV(a,b).  Also 

BCORR(a,b)  =  A(ao)  +  [A(bQ)  —  A(aQ-l-bQ)] 

ao=min{a,b},  bQ  =  max{a,b), 
where  (32)  and  (35)  are  applied. 

If  a<b,  then  BETALN(a,b)  can  be  accurately  computed  when  a>l.  If  a>8  then 
BETALN(a,b)  =  (.5ln  {‘In)  - .5  In  b)-(- BCORR(a,b)-u  - v 
11=  —  (a— .5)  ln[a/(a  +  b)^,  v  =  l)  ln(l  +a/b) 
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is  applied.  If  2<a<8  then  a  is  reduced  to  the  interval  [1,2]  by 
B(a,b)=^^^B(a-l,b). 

Consequently,  it  can  be  assumed  that  a<2  . 

If  b>8  then 

In  B(a,b)=ln  r(a)  +  ALG DI  V(a,b) 

is  applied.  If  2<b<8  then  b  is  also  reduced  to  the  interval  [1,2]  when  a>l.  Thus,  we  need  only 
consider  the  cases;  l<a<2,  1<  b<2,  or  a<l  and  b<8.  If  a>l  then 
In  B(a,b)=ln  r(a)  +  ln  r(b)-ln  r(a+b) 

is  appropriate.  No  loss  of  accuracy  due  to  subtraction  can  occur  since  Inr(a),  Inr(b),  and  — lnr(a  +  b) 
are  nonpositive.  However,  subtr2u:tion  does  occur  when  a<l.  It  currently  is  not  clear  how  loss  of 
accuracy  due  to  subtraction  can  be  avoided  when  a<l.  Therefore,  in  this  case  BETALN  is  not  used  in 
BRATIO. 

If  min{a,b}<8  then  BRCOMP(a,b,x,y)  can  be  computed  directly  from  its  definition. 
Otherwise, 

BRCOMP(a,b,x,y)  = 

z=^a</>(  1  —  A/a)  +  b0(  1  +  A/b)^  +  BCORR(a,b) 
is  used,  where  A  is  given  in  (15.1). 


ab  -z 
^25r(a  +  b) 
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V.  CONCLUDING  REMARKS 


Formulas  (12),  (13),  (14),  (15),  and  (16)  for  Ix(a,b)  are  of  the  form  7S,  where  S  is  a  series.  For 
example,  for  (15)  7  =  x^y^/B(a,b)  and  S  is  a  continued  fraction.  On  the  CDC  6000  —  7000  series 
computers,  almost  no  error  is  generated  in  computing  the  series  S  over  the  domains  specified  in  section 
III.  The  series  is  normally  accurate  to  within  1  or  2  units  of  the  14‘^  significant  digit.  However,  the 
precision  of  the  factor  7  is  restricted  by  the  inherent  error  of  lx(a,b).  Extensive  testing  on  the  CDC 
6000  —  7000  series  computers  comparing  the  results  obtained  by  BRATIO  with  results  from  double 
precision  code,  indicates  that  the  precisions  of  the  values  obtained  for  lx(a,b)  and  1— lx(a,b)  by 
BRATIO  approximate  the  inherent  errors  of  these  functions  up  to  a  maximum  of  14  significant  digits. 
On  any  computer,  accuracy  is  restricted  to  14  digits  because  of  the  algorithms  used  for  the  auxiliary 
functions  in  section  IV. 

On  the  CDC  6000  —  7000  series  computers,  a  maximum  of  7  terms  of  the  scries  (14)  for 
BGRAT,  and  a  maximum  of  11  terms  of  (16)  for  BASYM  were  observed  for  the  domains  specified  in 
section  III.  Frequently,  40  or  fewer  terms  of  (12)  for  BPSER  suffice,  but  a  maximum  of  92  terms  has 
been  observed  when  a  is  small,  b  is  large,  and  x  « .3.  Also,  40  or  fewer  terms  generally  suffice  for  the 
continued  fraction  (15),  BFRAC,  but  a  maximum  of  58  terms  has  been  observed  when  a  or  b  is 
exceedingly  large  and  x  2»a/(a  +  b). 

In  practice,  BRATIO  has  been  fouinl  to  l)e  a  reliable  and  efficient  subroutine.  As  was  noted  in 
the  previous  sections,  in  order  to  develop  such  a  subroutine  new  formulas  were  needed  for  ix(a,b),  a 
surprisingly  elaborate  specification  of  the  domains  of  usage  for  the  various  formulas  had  to  be  given, 
and  a  number  of  auxiliary  functions  had  to  be  treated  with  extreme  rare.  'I  hus,  the  development  of 
BRATIO  for  efficiently  computing  Ixfa.b)  and  lx(a,b)  to  high  relative  accuracy  was  not  a  simple 
task. 
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VI.  NUMERICAL  EXAMPLES 

A  collection  of  16  examples  using  DRATIO  is  given  in  Fig.2.  The  results  were  obtained  using 
the  CDC  6000-7000  series  single  precision  truncation  floating  arithmetic.  As  was  noted  in  section  III, 
the  function  BUP  is  used  only  with  BPSER  or  BGRAT  and  apjiears  in  (20),  (22),  (23),  and  (24).  The 
following  three  cases  illustrate  its  use.  The  quantity  c  below  is  set  to  approximately  .710E1-14  (see 

P-11). 

Case  (1).  a=.10,  b=14.5,  x=.29,  y=.71 
From  (20)  Wq  =BUP(14.5,  .10,  .71,  .29,  20,  e) 

=  1.7, (14.5,  .10)-I.7,(34.5,  .10) 

=  .17776  09989  0838E-3. 

Hence,  if  w  is  assigned  the  initial  value  wq  then  a  call  to  BGRAT(34.5,  .10,  .71,  .29,  w,  15f,  lERR) 
yields  the  value 

w  =Wo-f-1.7,(34.5,  .10)  =  I.7,(14.5,  .10) 

=  .17785  31648  7898E-3. 

Also 

I.29(.10,  14.5)=.99982  21468  3512. 

Case  (2).  a=1.5,  b  =  20.5,  x  =  .065,  y  =  .935 
(Note  that  bx>.70,  A  =  a  — (a-|-b)x  =  .07>0) 

From(22)  Wq  =BUP(.50,  1.5,  .935,  .065,  20,  e) 

=  I.935(.50,  1.5)-I.935(20.5,  1.5) 

=  .56745  07805  9439. 

Then 

1.065(1.5,  20.5)=Wo-l-BPSER(1.5,  .50,  .065) 

=  Wo-l-. 71754  32115  7741E-3 
=  .57462  62127  1016, 

1.935(20.5,  1.5>=.42537  37872  8984. 

Case  (3).  a=10.5,  b=1.5,  x  =  .80,  y  =  .20,  (A>0). 

From  (24)  Wq  =BUP(.50,  10.5,  .20,  .80,  I,  c) 

=  I.jo(-50,  10.5)-I.2o(1-5,  10.5) 

=  .15518  20005  6352, 

Wo  =  Wo-|-BUP(10.5,  .50,  .80,  .20,  20,  c) 

=  Wo  +  I.8o(10.5,  .50)-1.8o(30.5,  .50) 

=  Wo -(-.32149  19363  8971E-1 
=  .18733  11942  024.5. 
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Hence,  if  w  is  assigned  the  initial  value  Wq  then  a  call  to  BGRAT(30.5,  ,50,  .80,  .20,  w,  15f  ,IERR) 
yields  the  value 

w  =Wo  +  I.go(30-5t  -50) 
or 

I.8o(10.5,  1.5)=. 18756  94122  3880. 

Also, 

1.2o(1-5,  10.5)=.81243  05877  6120. 


a 

b 

X 

y 

1-Ix(a,b) 

.1 

.8 

.40 

.60 

.88776  70523  5302E+00 

.11223  29476  4698E+00 

.1 

.8 

.60 

.40 

.92957  83432  6833E+00 

.70421  65673  1668E-01 

.1 

2.3 

.40 

.60 

.97448  97683  7361E+00 

.25510  23162  6386E~01 

.1 

2.3 

.60 

.40 

.99196  58486  2884E+00 

.80341  51371  1598E-02 

5.0 

40.0 

.99 

.01 

.10000  00000  OOOOE+01 

.13053  04681  1410E-74 

5.0 

10.0 

.99 

.01 

.10000  00000  OOOOE+Ol 

.96509  74271  4997E-17 

10.0 

38.0 

.02 

.98 

.26944  43561  3309 E-07 

.99999  99730  5556E+00 

70.0 

10.0 

.85 

.15 

.23472  44941  6827E+00 

.76527  55058  3173E+00 

*  70.0 

50.0 

.99 

.01 

.10000  00000  OOOOE+Ol 

.54279  07073  1686E-66 

.  70.0 

50.0 

.10 

.90 

.47438  77486  2163E-38 

.10000  00000  OOOOE+Ol 

*  75.0 

50.0 

.10 

.90 

.61550  21193  1591E-42 

.10000  00000  OOOOE+Ol 

♦  .500.0 

501.0 

.50 

.40 

.99999  99999  3299E+00 

.67009  77013  4757E-10 

500.0 

501.0 

.40 

.60 

.10148  03038  4399E-09 

.99999  99998  9852E+00 

1000.0 

1001.0 

.49 

.51 

.19153  11043  9543E+00 

.80846  88956  0457E+00 

1001.0 

1000.0 

.49 

.51 

.17957  42144  6754 E+00 

.82042  57855  3246E+00 

(a=:5.0E+20, 

b  = 

5UE+3,  x=1.0,  y=1.0E 

-17 

Ix(a,b)  =  . 49811  93659  6617,  Iy(b,a)  =  .50188  06340  3383) 

Due  to  inherent 

error,  the  4  starred  examples  are  correct  to  within  1  unit  of  the  12lh 

significant  digit. 

All  other  cases  arc  correct  to  within  5  units  of  the  14lh  significant  digit. 

Fig.  2  16  Examples  of  BRATIO 
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BOUNDS  ON  Ix(a,b)  WHEN  x=a/(a+b),  MIN{a,b}>l 


The  purf)ose  of  this  appendix  is  to  show  that 

1/2  <  Ip(a,b)  <  1— e~*  ifl<a<b 

l/e<Ip(a,b)<  1/2  ifl<b<a 


(A-1) 


where  p=a/(a+b)  and  q=  1  — p=b/(a  +  b).  Since  I_5(a,a)  =  .5  from  (4),  we  shall  assume  that  a^b. 

Lemma  (A-1):  If  A(h)  =  (l -Hl/h)"^  for  h>0  then  A  is  a  decreasing  function.  Also,  A-»l  when 
h-*  0  and  A-»l/e  when  h-»  oo. 

The  lemma  is  given  for  reference.  It  is  a  well  known  result. 

Corollary  (A-1):  If  a  >  1  then  1/c  <  Ip(a,l)  <  1/2. 

Proof:  From  (1)  we  obtain  Ip(a,l)=p^  =  (l-)-l/a)’^.  Then  from  the  lemma  the  corollary 

follows. 


Corollary  (A-2):  If  b  >  1  then  1/2  <  Ip(l,b)  <  1  —  1/e. 

Proof:  Immediate  from  (4)  and  corollary  (A-1). 

Since  the  above  corollaries  hold,  in  order  to  verify  (A-1)  it  suffices  to  assume  that  a,  b>l.  The 


following  reasoning  is  due  to  James  C.  Perry  (NSWC).  Let 


B=|‘t®-*(l-t)^-Mt 

(A-2) 

Bp=|^t^-‘(l-t)^-Mt 

(A-3) 

Bp=|V-‘(l-t)^-‘dt. 

(A-4) 

Ifx=(l  — t)^^*  then 

(A-5) 

Bp'll' 

(A-6) 

where 

A=q'’''*  =  (l-|-a/b)~^'*. 

Also  1/e  <  A  <  1  by  lemma  (A-1),  and  we  now  consider  the  functions 

(A-7) 

f(x)  =  x-x‘^*^^  and  F(x)  =  f(x)*-‘ 

(A-8) 

for  0  <  X  <  1 . 
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From  (A-7)  and  (A-8)  we  note  that 

f(0)=F(0)=0,  f(l)  =  F(l)=0, 

f(x)  and  F(x)  have  a  unique  maximum  at  x  =  A,  and 

f  ^(x)  and  F^(x)  are  positive  (negative)  for  0  <  x  <  A  (A  <x  <  1). 

For  any  X2  such  that  A  <  Xg  <  1  let 

t(x,)  =  (A-9) 

where  Xj  is  the  unique  value  in  [0,A)  where  f(xj)=f(x2)  [see  Fig.  3]. 

We  now  examine  r(x2). 


Lemma  (A-2);  If  g(x)  is  twice  continuously  differentiable  on  [0,d],  g(0)=g^(0)  =  0,  and  g'*  is 
increasing  (decreasing)  on  (0,d],  then  h  is  continuous  and  increasing  (decreasing)  on  [0,d]  where 

(  g(x)/x^  (x>0) 

h(x)  =< 

I  g"(0)/2  (x=0). 

Proof:  Assume  that  g"  is  increasing  on  (0,d].  Since  h^(x)  =  (xg^  — 2g)/x^  for  x  >  0,  let 
k(x)=xg'  — 2g  for  x>0.  If  x>0  then  there  exists  0<^<x  such  that  g'(x)=g"(Ox  by  the  Mean- 
value  theorem,  so  that 

k'(x)  =xg'\x)-g'(x)  =  x[g''(x)-g"(^)]  >  0. 

Therefore,  k(x)  is  increasing  for  x  >  0.  Since  k(0)  =  0,  h'(x)=-k(x)/x^  >  0,  which  proves  that,  h  is 
increasing  on  (0,d].  Also,  since  g(0)  =  g'(0)  =  0,  h(x)-*g”(0)/2  when  x-»0  by  L’Hopital’s  rule.  If  g”(x) 
is  decreasing  on  (0,d]  then  apply  the  lemma  to  — g. 

Theorem  (A-1):  Let  r(x)  denote  the  function  defined  by  (A-9).  If  a<b  (a>b)  then  r(x)  is 
increasing  (decreasing)  for  A  <x  <  1.  Also,  r-»l  when  x-»A. 

Proof:  If  <6(h)  =  f(A)-f(A4-h)  for  0<h<l-A  then  (^(0)  =  0,  ^'(0)= -f'(A)  =  0,  and 

(^"(h)  =  (a/bq)(A-f-h)  for  h>0.  Hence,  if  a<b  (a>b)  then  <i>"  is  decreasing  (increasing),  so  that 


A-2 


(^(h)  =  (^(h)/h^  is  decreasing  (increasing)  by  lemma  (A-2).  Also,  (^(h)-»(^^^(0)/2=  a/(2bA)  when  h-»0 
b/a 

by  the  lemma  and  A=q 

Similarly,  if  V’(h)=f(A)— f(A  — h)  for  0<h<A  then  V’(h)  =  V’(h)/h^  is  increasing  (decreasing) 
when  a<b  (a>b),  and  ^(h)-nl’^^(0)/2=a/(2bA)  when  h-»0. 

Now  if  A  <x  <1  and  x,  <  A  where  f(x2)=f(x),  then 


,2_(x-AV_f(A)-f(xi)  (x-A)^  _t/-(A-X2) 
-U-Xi;  (A-x,)2  f(A)-f(x)~  ^(x-A) 


so  that  “TT^r^  =  1  when  x-»A.  If  a  <  b  then  when  x  increases.  A— Xj  increases,  i/i(A— Xj)  increases, 

3./  y£UA  f 

S(x  — A)  decreases,  and  hence  r^  increases.  Similarly,  r^  decreases  when  a>b.  Consequently,  since  r>0 
the  theorem  follows. 


Theorem  (A-2):  If  1  <  a  <  b  then  Ip(a,b)  <  1  —  1/e. 

Proof:  Since  a  <  b,  r(x)  <  r(l)=  (1  —  A)/A  for  A<x<l  by  theorem  (A-1).  Also,  since 
Ip(l,b)  <  1  —  1/e  by  corollary  (A-2),  we  shall  assume  that  a>  1.  Now  let 

A<x<l  (A-10) 

so  that  0  <  t  <  A.  For  x  >  A  let  Xi  <  A  where  f(xi)=  f(x).  Then  t  <  A,  and  from 

r(x)  =  v^~-^  <  we  have  x,  <  t.  Therefore, 

A  — Xi  —  A  ‘  ” 

F(l-l^t)=F(x)=F(x2)<F(t) 

so  that 

Bp=|j'F(x)dx  =  |4ij^^F(l-4il)dt 
^5  ^lo  (i-l)Bp<(e-l)Bp 


from  (A-5)  and  (A-6).  Therefore,  from  (A-2)  — (A-4),  Bp  <  (e— l)Bp  =(e— 1)(B  — Bp)  so  that 
lp(a,b)  =  Bp/B<  1-1/e. 

Theorem  (A-3);  If  1  <  b  <  a  then  Ip(a,b)  <  1/2. 

Proof:  Since  a  >  b,  r(x)  <  1  for  A  <x  <  1  by  theorem  (A-1).  Also,  since  Ip(a,l)  <  1/2  by 
corollary  (A-1),  we  shall  assume  that  b  >  1.  Now  let 

t=2A-x  A<x<l.  (A-11) 

Since  b/a  <1,  .5  <  A  <  1  from  (A-7)  and  lemma  (A-1).  Therefore,  0  <  2A  —  1  <  A  and  2A  —  1  <  t  <  A. 
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For  x>A  let  Xj  <  A  where  f(xi)=f(x).  Then  t  <  A,  and  from  t(x)=y — ^<1  we  have  Xj  <  t. 
Therefore, 

F(2A-t)=  F(x)=  F(xi)  <  F(t) 

so  that 

from  (A-5)  and  (A-6).  Therefore,  from  (A-2)— (A-4)  Bp  <  Bp=  B  — Bp  so  that  Ip(a,b)=  Bp/B  <  1/2. 

Theorem  (A-4):  If  a,b  >  1  then 

1/2 < Ip(a,b)  <  1  —  1/e  whena<b 
l/e<  Ip(a,b)  <  1/2  whena>b. 

Proof:  Immediate  from  (4)  and  theorems  (A-2)  and  (A-3). 


A-4 


APPENDIX  B 


MQNQTONIC  PROPERTIES  OF  Ix(a,b) 


MONOTONIC  PROPERTIES  OF  Ix(a,b) 


The  purpose  of  this  appendix  is  to  show  that  Ij((a,b)  is  a  decreasing  function  of  a  and  an 
increasing  function  of  b. 

Theorem  B:  If  a>0,  b>0,  and  0<x<l,  then 

~  ~  (7a  — 

Proof:  Let 

f(t;a,b)  =  t^“‘(l-t)^“', 
hereafter  simply  denoted  by  f(t).  We  have 
Ix(a,b)  =  Bx(a,b)/B(a,b) 


Then 


Bx(a,b) 


’'t'‘-‘(l-t)^  ‘dt, 
0 


(i-t)^''dt. 


aix(a,b) 

’‘u^-'(l-u)‘’"'ln  udu-Bx(a,b) 

0 

't^-'(l-t)^-'lntdt 

0 

5a  “ 

B(a,b)  B(a,b) 

[B(a,b)J  =  ['  r 1  -t)'>-'  ua-'(l  -u)*^-’  (In  u ~ln  t)du 


I  0 
rx  rx 


dt 


'X  rx  ri  rx 

=  [^f(t)f(u)()i)ij— )))t)dujdt+  [^f(t)f(u)  (In  u  — In  t)dujdt. 


But  the  first  double  integral  on  the  right  is  zero.  Indeed,  with 


•X 

0 


[^f(t)f(u)(ln  n— In  t)du^dt. 


where  the  subscripts  are  u.sed  to  indicate  the  order  of  integration  (u  then  t),  we  have,  by  renaming  the 
variables, 

u,t  ‘t,u' 

Interchanging  the  order  of  integration  then  gives 

T  — 'P  —  —  T 

‘u,l~‘t,u“  ‘u,t’ 

which  implies  that  T^  ^=0-  'I'herefore, 

[B(a,b)7  ^^^^^=|^|^(t)f(u)(lnu-lnt)dudt,  (H-l) 

which  is  negative  since  Inu  — lnt<()  for  u<t.  Thus, 

— for  all  a>(),  b>0,  0<x<l. 

/la  —  —  — 
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Remarks.  A  more  direct  proof  of  (B-1)  can  be  obtained  by  differentiating  1/[1+Iy(b,a)/Ix(a,b)].  Also, 
since  In u  — In t<u  — 1<0  for  0<u<t<l,  from  (B-1)  we  obtain  the  slightly  stronger  inequality 

-\l\l 

B(a,b)Bx(a+l,b)  B(a+l,b)Bx(a,b) 

[B(a,b)]2  [B(a,b)]2 

r+b  ^  ~ 

Corollary:  Ifaj<a2  then  Ix(a2,b)<Ix(a|,b). 

Corollary:  If  bj<b2  then  Ix(a,bi)<Ix(a,b2). 

Proof:  Since  Iy(b,a)  decreases  as  b  increases,  lx(a,b)  =  1  — Iy(b,a)  increases  as  b  increases. 

We  close  by  noting  that  precisely  the  same  technique  used  in  Theorem  B  establishes  a 
corresponding  result  for  the  incomplete  gamma  function.  If 

rx 

P(a,x)=pf-^  e-H®-'dt,  a>0,  x>0, 

then  and  hence  if  ai<a2then  P(a2,x)<P(ai,x). 
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UNIFORM  ASYMPTOTIC  BEHAVIOR  OF  (16) 


The  purpose  of  this  appendix  is  to  show  that  (16)  has  uniform  asymptotic  behavior  for  a<b 
and  a(l+a/b)  oo.  The  proof  will  depend  on  showing  that  the  n**  derivative  of  the  function  f(v)=:v/s 
(v>0)  hcis  a  bound  not  depending  on  h=a/b  when  a<b  and  n  is  odd. 

Theorem  (C-1);  If  a<b  then  0<  f^(v)<l  for  v>0  . 

Proof:  Since  ^=|(1— s)(l +hs)  and  f=g  we  have 


Also 


ds_l_yf 

dv  ®  §3 


(1— s)(l+hs). 


^  (l-s)(H-hs)=  s" 

s  n>i 


An_  1+h"^*  h+h" 
2q  (n  +  l)(n  +  2)^  n(n+l) 


n  even 


An  _  1-h"'*'*  I  h-h"  yn 

2q  (n  +  l)(n+2)^  n(n  +  l)  ~ 


n  odd. 


(C-1) 


(C-2) 


(C-3) 


Therefore  f^(v)  =  52  An  s"  ‘>0  is  an  increasing  function  for  s>0  (and  hence  for  v>0).  From  (C-2) 
n>i 

1  — 53Ans">0  for  0<s<l,  so  that  X)  ^n^l-  Consequently,  f ^(v)  =  ^  An  s”“' <  5lAn<l  for 
n>i  n>i  n>i  n>i 

v>0. 

Theorem  (C-2);  If  a<b  then  f^(v)  -»  I  when  v  -»  oo. 

Proof:  From  (C-3)  An/(2q)  =  (l-l-h)  ^  l/[n(n-t- 1)],  where 

n>i  n>2 

l/[n(n-M)]=  [l/n  -l/(n-l-l)]=  1. 
n>i  n>i 

From  (16.7)  we  have 

v2=jp^(-ln(l-s)-i|n  (l-l-hs)].  (C-4) 

Consequently,  we  now  consider  v  as  a  function  defined  for  0<h<l  and  0<s<l. 

Lemma  (C-1):  If  n=0,l,2,...  then  v"(l— s)  can  be  extended  to  a  continuous 
function  ^n(li<^)  fo*’  0<h<l  and  0<s<l.  Also,  ^n(h>l)=0  for  0<h<l. 

Proof:  From  (C-i)  we  obtain 
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Ojj(s)=[-ln  (l-s)]*'  (l-s)^  0<s<l 

g(z)=  — In  (1 +z)/z  0<z<l. 

By  L’Hopital’s  rule  g(z)-*  — 1  when  z-»0.  Hence,  if  g(z)=g(z)  for  0<z<l  and  g(0)  =  — 1,  then  g(z)  is 
continuous  for  0<z<l.  Also  if  w=—  In  (1— s)  then 
lim  a,(s)=lim  w*'exp(-2w)  =  0. 

Thus,  if  Qj^(s)  =  aj^(s)  for  0<s<l  and  5|^(1)  =  0,  then 

l/2 

is  a  continuous  extension  of  v"(l  — s)  for  0<h<  1  and  0<s<  1.  Also,  (^f)(hT)  =  0  since  each  5j^(l)  =  0. 

Theorem  (C-3):  If  n>2  then  f^"\v)  can  be  extended  to  a  function  that  is  continuous  for 

0<h<l  and  0<8<1.  Also,  f^”\v)  -♦  0  when  v  -*  oo. 

Proof:  Given  a  fixed  value  0<h<l.  Then  for  any  polynomial  P  in  s  (and  h),  k  =  0,l,2,...  , 
and  m  =  l,2,...  , 

^  [34(l-s)Pl=k  (1-s)  (l-s)^(l+hs)P 

-  4^  (l-s)(l+hs)P  +  4^  (l-s)2(l+hs)P' 
s  s 

where  P^  is  the  derivative  dP/ds.  By  induction  (starting  with  C-1)  it  follows  for  n>2  that  f*^*'\v)  is  a 
finite  sum  of  terms  of  the  form  (v^/s*’^)(l  — s) P.  Consequently,  by  Lemma  (C-1)  r^"\v)  has  a  con¬ 
tinuous  extension  V'n(h,s)  for  0<h<l  0<s<l.  Also,  V’n(h-1)  =  0  so  that  f*'''(v)  -»  0  when  v-*cc. 

Remark.  For  n>2,  since  f(v)  =  53 ^  neighborhood  of  0  and  f'"'(v)  -♦  0  when  v-*oc,  by 

k>o 

Theorem  (C-3)  f*”*(v)  is  bounded  for  v>0.  However,  the  bound  may  depend  on  h. 

We  now  prove  that  a  bound  exists  not  depending  on  h  when  n  is  odd. 

Definition:  Given  a  series  53  where  each  On(l')  continuous  for  0<h<l. 

n>o 

Let  Mn=max{  |ajj(h)|  :0<h<l  }.  'I'lien  the  series  53onS*'  is  said  to  satisfy  condition  (»)  when 

53  M^s  converges  for  |s|  <  1 . 
n>o 
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Lemma  (C-2):  f'(v)  =  ^  *  satisfies  condition  (♦). 

n>i 

Proof:  From  (C-3)  is  continuous  for  0<h<l  and 

I  A  |<; _ 4 _ I  4 

'  "'-(n  +  l)(n+2)^n(n  +  l)  • 

Also  l/[n(n  +  l)]=:l. 


f  k  ^ 

Theorem  (C-4):  For  k  odd,  the  series  Y!,  satisfies  condition  (*). 

n>o 

Proof:  By  Lemma  (C-2),  the  theorem  is  true  for  k  =  l.  Assume  inductively  that  it  is  true  for 

odd  k.  If  e[i=max{  lan(h)l  :  0<h<l  }  then  let  g(s)=  Y  Now 

n>o 

f'‘‘^'V)  =  (  Enons"-')|(l-s)(l-^hs)=|  ^/?n(h)s"  (C-5) 

^  n>i  '  n>o 

0Q  =  au  /^n=(n  +  l)an  +  !~n(l-h)on-h(n-l)ftn-i  (n>l)- 


Since  Q[i(h)  is  continuous  for  0<h<l,  therefore  0n(h)  is  continuous  for  0<h<l.  Also  |/?nl£^^n 

where  Mo=ei  and  Mn  =(n  + 1  )e(,  +  , -fnen -|-(n  — l)en_;  (n>l).  Since  Y  Mn®"  's  the  series 

for  g^(s)(l -l-s  +  s^),  Y  J^n®"  converges  for  |s|<l  and  Y  ^n  s”  satisfies  condition  (*). 
n>o  n>o 

From  (C-5)  we  obtain 

f'‘‘^  =  '(v)=f'(v)  E  /?ns"+  En/?ns""'  i^'(l-s)(l +hs) 
n>o  n>i  s 

=  E  ^ns""  E  /?n«''  +  En/3ns"-‘  (l-E  ^ns"  ). 

n>i  n>o  n>i  '  n>i  ' 

Consequently,  f**^^^*(v)  =  5!^  ifi„(h)s’^  where 

n>o 

6o=0\  +  00^1 

^n=(n  +  l)/fn  +  i  +  /^o^n*i-'E  A;._i  (">!)•  (C-G) 

i=o 

Then  is  continuous  for  0<h<  1.  Also,  if  c,,  =  niax  (  A,,  :  0<h<  1  }  then 


I'^nlS  ('>+*)  ^*11  +  1  +  M()Cn  +  i  +  Np 

Nn  =  E'iMi^,c^^_.  (n>l),  N„  =  0. 

i=o 

Now  E  (o  O'^^n  +  i  •'>"  converges  for  |s|<l  since  it  is  the  derivative  of  E^'ns'*> 
n>o  n>o 

and  E*^n  +  l^^  converges  for  |s|<l  by  Lemma  ('-2.  Also 
n>o 


C-3 


E  Nnlsl"  =  E  nM„„|sr  E  Cnhl"  <  E  (n  +  l)Mn,ilsl"  E  Cn|s|", 

n>o  n>o  n>o  n>o  n>o 

which  converges  for  lsl<l.  Consequently,  =  E  ^n®”  satisfies  condition  (+). 

n>o 

(k) 

Theorem  (C-5):  If  k  is  odd  then  f  (v)  can  be  extended  to  a  function  that  is  continuous  for 

ik) 

0<h<l  and  0<s<l.  Consequently,  |f^  ^(v)!  is  bounded  for  v>0  by  a  value  that  does  not  depend  on  h. 


Proof:  It  has  already  been  verified  that  f  \v)  can  be  extended  to  a  function  that  is 

f  k ) 

continuous  for  0<h<l  and  0<s<l.  Let  0<So<l.  Since  f^  (v)  =  E  ^n(h)s''  satisfies  condition  (*), 

n>o 

E  converges  for  M[i=rnax{  |ajj(h)|  :  0<h<l}.  Therefore,  by  the  Weierstrass  M-test, 

n>o 


E  ans”  converges  uniformily  for  0<h<l  and  0<s<So.  Since  each  On(h)s"  is  continuous  for 
n>o 

0<h<l  and  0<s<So,  therefore  E  ®n  continuous  for  0<h<l  and  0<s<Sq. 

n>o 
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FORTRAN  LISTING  OF  BRATIO 


In  this  appendix  the  code  for  BRATIO  is  given.  The  following  functions  arc  used. 
ERF(x)=erf  x 


,  erfc  X 

ERFCl(i,x)  =  J 

^  exp(x^)erfc  x 

REXP(x)  =  cxp(x)-l 
ALNREL(a)  =  ln(l+a) 
RLOGl(a)=a-ln(l+a) 

GAMl(a)  =  l/r(l+a)-l 
GAMLN(x)  =  ln  r(x) 

GAMLNl(a)  =  ln  r(l+a) 
ALGDIV(a,b)  =  ln[r(b)/r(a  +  b)] 
BCORR(a,b)  =  A(a)  +  A(b)-A(a  +  b) 
BETALN(a,b)  =  ln  B(a,b) 
BRCOMP(a,b,x,y)=x^y^/B(a,b) 


if  i=0 
if  i^O 

a>  —  1 
a>  —  1 
—  .5<a<  1 .5 
x>0 

-.2<a<  1.2.5 
a>0,  b>8 

a,b>8 
a,b>0 

a,b>0,  0<x<l,  y=l— X 


'I'hese  functions,  written  by  A. 11. Morris,  are  part  of  the  NSWC  mathematics  subroutine  library  [9]. 
Machine  -  Dependent  Constants 

'I'he  function  SPMPAR  provides  the  machine-dependent  constants  needed  by  BRATIO.  It  is 
ncces.sary  that  SPMPAR  be  properly  (kTincd  for  the  computer  arithmetic  being  used.  The  constants  are 
defined  in  the  in-line  documentation  of  SPMAR.  Values  for  these  constants  are  given  in  the  in-line 
documentation.  Sl’MPAR,  released  by  Argonne  National  Laboratory,  is  an  ada[)tation  of  the  Bell 
Laboratories  function  III  MACH  [Dl]. 

Transportability 

All  coding  adheres  to  the  1966  and  1977  ANSI  FORTRAN  standards.  It  is  a.ssumed  that  a 
floating  point  arithmetic  of  6  or  more  digits  is  being  used.  I'he  codes  were  designed  specifically  for 
k-digit  arithmetics  where  If  k>l4  then  only  li-digil  accuracy  will  normally  be  obtained. 

Reference 

Dl.  Fox,  P.A.,  Hall,  A.D.,  and  Schryer,  N.L.  The  PORT  Mathematical  Subroutine  Library, 
ACM  Trans.  Math  Software  4  (1978),  pp.  101-126. 
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SUBROUTINE  BRATIOIA.  B,  X,  Y,  W.  W1.  lERR) 

C - - - - - 

c 

C  EVALUATION  OF  THE  INCOMPLETE  BETA  FUNCTION  IX(A,B) 

C 

C  - - 

C 

C  IT  IS  ASSUMED  THAT  A  AND  B  ARE  NONNEGATIVE,  AND  THAT  X  . LE .  1 

C  AND  Y  =  1  -  X.  BRATIO  ASSIGNS  W  AND  W1  THE  VALUES 

C 

C  W  =  IX(A,B) 

C  W1  =  1  -  IX(A,B) 

C 

C  lERR  IS  A  VARIABLE  THAT  REPORTS  THE  STATUS  OF  THE  RESULTS. 

C  IF  NO  INPUT  ERRORS  ARE  DETECTED  THEN  lERR  IS  SET  TO  0  AND 

C  W  AND  W1  ARE  COMPUTED.  OTHERWISE.  IF  AN  ERROR  IS  DETECTED, 

C  THEN  W  AND  W1  ARE  ASSIGNED  THE  VALUE  0  AND  lERR  IS  SET  TO 

C  ONE  OF  THE  FOLLOWING  VALUES  . . . 

C 


c 

lERR 

1 

IF 

A 

OR  B 

IS  NEGATIVE 

c 

I  ERR 

= 

2 

IF 

A 

=  E 

=  0 

c 

lERR 

a 

3 

IF 

X 

.  LT  . 

0  OR  X  . GT . 

c 

lERR 

= 

4 

IF 

Y 

.  LT  . 

0  OR  Y  .GT  . 

c 

lERR 

= 

5 

IF 

X 

+  Y 

.NE.  1 

c 

lERR 

= 

6 

IF 

X 

=  A 

=  0 

c 

lERR 

a 

7 

IF 

Y 

=  B 

=  0 

C 

C . 

C  WRITTEN  BY  ALFRED  H.  MORRIS.  JR. 

C  NAVAL  SURFACE  WEAPONS  CENTER 

C  DAHLGREN.  VIRGINIA 

C  REVISED  . . .  JUNE  1988 

C . 

REAL  LAMBDA 


EPS  IS  A  machine  DEPENDENT  CONST  ANT .• E PS  IS  THE  SMALLEST 
FLOATING  POINT  NUMBER  FOR  WHICH  1.0  +  EPS  .GT.  1.0 


C 


C 


C 


EPS  =  SPMPARI  1  ) 


W  =  0.0 
W1  =  0.0 

IF  (A  .LT.  0.0  .OR.  B  .LT.  0.0)  GO  TO  300 
IF  (A  .EO.  0.0  .AND.  B  .EO.  0.01  GO  TO  310 
IF  (X  .LT.  0.0  .OR.  X  .GT.  1.0)  GO  TO  320 
IF  (Y  .LT.  0.0  .OR.  Y  . GT .  1.0)  GO  TO  330 
Z  =  DBLE(X)  +  DBLE(Y)  -  1 . DO 

IF  (ABS(Z)  .GT.  EPS)  GO  TO  340 


lERR  =  0 
IF  (X  .EO. 
IF  (Y  .EO. 
IF  (A  .EO. 
IF  (B  .EO. 


0.0)  GO  TO  200 
0.0)  GO  TO  210 
0.0)  GO  TO  211 
0.0)  GO  TO  201 


IND  =  0 
AO  =  A 
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c 

c 

c 


c 


c 


c 

c 

c 


c 


c 

c 

c 


BO  =  B 
XO  =  X 
VO  =  Y 

EPS  =  AMAXKEPS.  1  .E-15) 

IF  (AMINKAO,  BO)  .  GT .  1.0)  GO  TO  30 

PROCEDURE  FOR  AO  -LE.  1  OR  BO  .LE.  1 

IF  (X  .LE.  0.5)  GO  TO  10 

IND  =  1 

AO  =  B 

BO  =  A 

XO  =  Y 

YO  =  X 

10  IF  (AMAXKAO.  BO)  .GT.  1.0)  GO  TO  20 

IF  (AO  .GE.  AMIN1(0.2.  BO))  GO  TO  100 

IF  (XO*»AO  .LE.  0.9)  GO  TO  100  t 

IF  (XO  .GE.  0.3)  GO  TO  110 
N  =  20 
GO  TO  130 

20  IF  (BO  .LE.  1.0)  GO  TO  100 

IF  (XO  .GE.  0,3)  GO  TO  110 

IF  ( XO  .GE .  0.1)  GO  TO  21 

IF  ((XO*BO)»*AO  .LE.  0.7)  GO  TO  1(X) 

21  IF  (BO  .GT.  15.0)  GO  TO  131 
N  »  20 

GO  TO  130 

PROCEDURE  FOR  AO  .GT.  1  AND  BO  . GT .  i 

30  IF  (A  .GT.  B)  GO  TO  31 

LAMBOA  =  A  -  (A  +  B)*X 
GO  TO  32 

31  LAMBDA  »  (A  +  E)*Y  -  B 

32  IF  (lambda  GE.  0.0)  GO  TO  40 
IND  =  1 

AO  =  E 
BO  =  A 
XO  =  Y 
YO  =  X 

LAMBOA  =  A6S( LAMBOA) 

40  IF  (BO  .LT.  40.0  .AND.  BO*XO  .LE.  0.7)  GO  TO  100 
IF  (BO  .LT.  40.0)  GO  TO  140 
IF  (AO  .GT.  BO)  GO  TO  50 

IF  (AO  LE.  100.0)  GO  TO  120 
IF  (lambda  .GT.  0.03*A0)  GO  TO  120 
GO  TO  180 

50  IF  (BO  .LE.  100.0)  GO  TO  120 

IF  (LAMBOA  .GT,  0.03*80)  GO  TO  120 
GO  TO  180 

EVALUATION  OF  THE  APPROPRIATE  ALGORITHM 
100  W  =  BPSERIAO.  BO.  XO.  EPS) 
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c 


c 


c 


c 


c 


c 


c 

c 

c 


c 


c 


W1  =  0.5  +  (0.5  -  W) 

GO  TO  220 

110  W1  =  BPSER(BO.  AO.  YO,  EPS) 

W=0.5+(05-W1) 

GO  TO  220 

120  W  =  BFRAC(A0,  BO,  XO.  YC.  LAMBDA,  15.0*EPS) 

W1  =  0.5  +  (0.5  -  W) 

GO  TO  220 

130  W1  =  BUP(BO.  AO,  YC.  XO,  N.  Ef=S) 

60  =  BO  +  N 

131  CALL  BGRAT(BO,  AO,  YO ,  XO.  W1,  15.0*£PS,  IERR1) 

W  =  0.5  +  (0.5  -  W1  ) 

GO  TO  220 

140  N  =  BO 

BO  =  BO  -  N 

IF  (BO  .NF.  0.0)  GO  TO  141 
N  -  N  -  1 
BO  =  1.0 

141  W  =  BUP'SO,  AO.  YO.  XO.  N.  EPS) 

IF  (XO  .GT.  0.7)  GO  TO  150 

W  *  W  +  BPSER(A0,  BO.  XO.  EPS) 

W1  =  0.5  +  (0.5  -  W) 

GO  TO  220 

150  IF  (AO  .GT.  15.C)  GO  TO  151 

N  =  20 

W  =  W  +  BUP(AO.  BO.  XO.  YO.  N.  EPS) 

AO  =  AO  +  N 

151  CALL  BGRAKAO,  BO.  XO.  YO.  W.  15.0»£PS.  IERR1) 
W1  *  0.5  +  (0.5  -  W) 

GO  TO  220 

180  W  =  BASYM(AO.  BO,  XO,  YO.  LAMBDA.  100.0*EPS) 

W1  -  0.5  +  (0.5  -  W) 

GO  TO  220 

TERMINATION  OF  THC  PROCEDURE 

200  IF  (A  .EO.  0.0)  GO  TO  350 

201  W  =  0.0 
W1  =  1.0 
RETURN 

210  IF  (B  .EO.  0.0)  GO  TO  360 

211  W  =  1.0 
W1  =  0.0 
RETURN 

220  IF  (IND  .EO.  O)  RETURN 
T  =  W 
W  =  W1 
W1  =  T 
RETURN 
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c  ERROR  RETURN 

C 

300  I  ERR  =  1 
RETURN 

310  lERR  =  2 
RETURN 

320  lERR  =  3 
RETURN 

330  lERR  =  4 
RETURN 

340  I  ERR  '  5 
RETURN 

350  lERR  =  6 
RETURN 

360  I  ERR  -  7 
RETURN 
END 
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REAL  FUNCTION  BPSERlA,  B,  X.  EPS) 


POWER  SERIES  EXPANSION  FOR  EVALUATING  IX(A.B)  WHEN  B  . LE .  1 

OR  B»X  .LE.  0.7.  EPS  IS  THE  TOLERANCE  USED. 


REAL  N 

BPSER  =  0.0 

IF  (X  .EO.  0.0)  RETURN 


COMPUTE  THE  FACTOR  X* •A/ ( A»BETA ( A , B ) ) 


AO  =  AMINI(A.B) 

IF  (AO  .LT.  1.0)  GO  TO  10 

Z  =  A*ALOG(X)  -  BETALN(A.B) 

BPSER  =  EXP(Z)/A 
GO  TO  70 

10  BO  =  AMAXKA.B) 

IF  (BO  .GE.  8.0)  GO  TO  60 
IF  (BO  .GT.  1.0)  GO  TO  40 

PROCEDURE  FOR  AO  . LT .  1  AND  BO  .LE.  1 

BPSER  =  X**A 

IF  (BPSER  .EO.  0.0)  RETURN 
APB  =  A  +  B 

IF  (APB  .GT.  1 .0)  GO  TO  20 
Z  =  1.0+  GAMI(APB) 

GO  TO  30 

20  U  *  DBLE(A)  +  OBLE(B)  -  1 . DO 
Z  *  (1.0+  GAM1 (U ) )/APB 

30  C  =  (1.0  +  GAM1(A) )•( 1 .0  +  GAM1(B))/Z 
BPSER  *  BPSER»C*(B/APB) 

GO  TO  70 

PROCEDURE  FOR  AO  . LT .  1  AND  1  . LT .  BO  . LT .  8 

40  U  =  GAMLN1 ( AO) 

M  =  BO  -  1.0 

IF  (M  . LT .  1 )  GO  TO  50 

C  =  1 .0 
DO  41  I  =  1  ,M 
BO  =  BO  -  1.0 

41  C  =  C*(B0/( AO  +  BO) ) 

U  =  ALOG(C)  +  U 

50  Z  =  A*ALOG(X)  -  U 
BO  =  BO  -  1.0 
APB  =  AO  +  BO 

IF  (APB  .GT.  1 .0)  GO  TO  51 
T  =  1.0+  GAM 1( APB) 

GO  TO  52 

51  U  =  OBLE(AO)  +  DBLE(BO)  -  1 . DO 
T  =  (1.0+  GAMKU)  )/APB 

52  BPSER  =  EXP(Z)*(AO/A)*( 1.0  +  GAM1(B0))/T 
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GO  TO  70 

PROCEDURE  FOR  AO  . LT .  1  AND  BO  . GE .  8 

60  U  =  GAMLNKAO)  +  ALGO  I  V(  AO .  BO ) 

2  =  A»ALOG(X)  -  U 
BPSER  =  (A0/A)»EXP(2) 

70  IF  (BPSER  .EO.  0.0  .OR.  A  .LE.  0.1’EPS)  RETURN 


COMPUTE  THE  SERIES 


SUM  =  0.0 
N  =  0.0 
C  =  1.0 
TOL  =  EPS/A 
100  N  =  N  +  1 . 0 

C  =  C-(0.5  +  (0.5  -  B/N)  )»X 
W  =  C/(A  +  N) 

SUN.  =  SUM  +  W  . 

IF  (Al3S(W)  .GT.  TOL)  GO  TO  100 
BPSER  =  BPSER»(1.0  +  A«SUM) 

RETURN 

END 
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REAL  FUNCTION  BUP(A,  B.  X.  Y.  N,  EPS) 


EVALUATION  OF  IX(A,B)  -  IX(A+N.B)  WHERE  N  IS  A  POSITIVE  INTEGER. 
EPS  IS  THE  TOLERANCE  USED. 


REAL  L 

BUP  =  BRCOMP(A.B.X. Y)/A 

IF  (N  .EO.  1  .OR.  BUP  .EO.  0.0)  RETURN 

NM1  =  N  -  1 

APB  =  A  +  B 

AP1  =  A  +  1.0 

SUM  =  1.0 

D  =  1.0 

LET  K  BE  THE  INDEX  OF  THE  MAXIMUM  TERM 


K  =  0 

IF  (B  . LE .  1.0)  GO  TO  30 
IF  ( Y  .GT.  1  . E-4 )  GO  TO  10 
K  =  NM1 
GO  TO  20 

10  R  =  (B  -  1 .0)»X/Y  -  A 

IF  (R  .LT.  1.0)  GO  TO  30 
K  =  NM1 

IF  (R  . LT.  FL0AT(NM1  )  )  K  =  R 

AOD  THE  INCREASING  TERMS  OF  THE  SERIES 

20  00  21  I  =  1 ,K 
L  =  I  -  1 

D  =  ((APB  +  L)/(AP1  +  L))*X«D 

21  SUM  »  SUM  +  0 

IF  (K  .  EO.  NMD  GO  TO  40 

AOD  THE  REMAINING  TERMS  OF  THE  SERIES 

30  KP1  =  K  +  1 

DO  31  I  =  KP1 .NM1 

L  =  I  -  1 

D  =  ((APB  +  L)/(AP1  +  L))*X*D 
SUM  =  SUM  +  0 

IF  (D  .LE.  EPS*SUM)  GO  TO  40 

31  CONTINUE 


TERMINATE  THE  PROCEDURE 

40  BUP  =  BUP»SUM 
RETURN 
END 
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REAL  FUNCTION  BFRAC(A,  B,  X,  Y,  LAMBDA,  EPS) 


CONTINUED  FRACTION  EXPANSION  FOR  IX(A.B)  WHEN  A,B  . GT .  1. 
IT  IS  ASSUMED  THAT  LAMBDA  =  (A  +  B)»Y  -  B. 


REAL  LAMBDA,  N 


BFRAC  =  BRCOMP(A,B,X,Y) 

IF  (BFRAC  .EO.  0.0)  RETURN 

C  =  1.0+  LAMBDA 
CO  =  B/A 

C1  =  1.0+  1 .O/A 
YP1  =  Y  +  1.0 

N  =  0.0 
P  =  1.0 
S  =  A  +  1 .0 
AN  =  0.0 
BN  =  1.0 
ANP1  =  1.0 
BNP1  =  C/Cl 
R  =  Cl/C 

CONTINUED  FRACTION  CALCULATION 

10  N  =  N  +  1 .0 
T  »  N  J 

W  *  N«(B  -  N)*X 
E  *  A/S 

ALPHA  »  (P*(P  +  C0)*£»E)*(W*X) 

E  =  (  1 .0  +  T)/(C1  +  T  +  T) 

BETA  =  N  +  W/S  +  E*(C  +  N+YP1) 

P  *  1.0+  T 
S  =  S  +  2.0 

update  AN.  BN,  ANP 1 ,  AND  BNP  1 

T  =  ALPHA*AN  +  BETA+ANP1 
AN  =  ANP1 
ANP1  =  T 

T  =  ALPHA+BN  +  BETA»BNP1 
BN  =  BNP1 
BNP1  =  T 

RO  =  R 

R  =  ANP 1 /BNP  1 

IF  (ABS(R  -  RO)  .LE.  EPS*R)  GO  TO  20 

RESCALE  AN.  BN,  ANPI,  AND  BNP1 

AN  =  AN/BNP  1 
BN  =  BN/BNP  1 
ANPI  *  R 
BNP1  =  1.0 
GO  TO  10 


0-  10 


c 

c 


TERMINATION 


20  BFRAC  =  BFRAC»R 
RETURN 
END 
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REAL  FUNCTION  BRCOMP(A,  B.  X.  V) 
EVALUATION  OF  X» • A • Y • •B/BET A ( A , B  ) 


REAL  LAMBDA,  LNX .  LNY 


CONST  =  1/S0RT(2*PI ) 


DATA  CONST/ . 398942280401433/ 

AO  =  AMINKA.B) 

IF  (AO  .GE.  8.0)  GO  TO  100 

IF  (X  .GT.  0.375)  GO  TO  10 
LNX  =  ALOG(X) 

LNY  =  ALNREL(-X) 

GO  TO  20 

10  IF  (Y  .GT.  0.375)  GO  TO  11 

LNX  =  ALNREL(-Y) 

LNY  =  ALOG(Y) 

GO  TO  20 

11  LNX  =  ALOG(X) 

LNY  =  ALOG(Y) 

20  IF  (AO  .LT.  1.0)  GO  TO  30 

2  =  (A*LNX  +  B*LNY)  -  BETALN(A.B) 
BRCOMP  =  EXP(Z) 

RETURN 


PROCEDURE  FOR  A  . LT .  1  OR  B  . LT .  1 


30  BO  =  AMAXKA.B) 

IF  (BO  .GE.  8.0)  GO  TO  80 
IF  (BO  .GT.  1 .0)  GO  TO  60 

algorithm  FOR  BO  .LE.  1 

BRCOMP  =  EXP(A-LNX  +  B»LNY) 

IF  (BRCOMP  .EO.  0.0)  RETURN 

APB  =  A  +  B 

IF  (APB  .GT.  1 .0)  GO  TO  40 
2  =  1.0+  GAM1 ( APB ) 

GO  TO  50 

40  U  =  OBLE(A)  +  DBLE(B)  -  1.00 

2  =  (1,0  +  GAMKU)  )/APB 

50  C  =  (1.0  +  GAMKA  )  )•(  1  .0  +  GAM1(B))/2 
BRCOMP  =  BRCOMP* ( AO*C )/( 1 .0  +  AO/BO) 
RETURN 


ALGORITHM  FOR  1  . LT .  80  .LT.  8 

60  U  =  GAMLNKAO) 

N  =  BO  -  1.0 

IF  (N  , LT .  1  )  GO  TO  70 

C  =  1.0 
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61 


DO  61  I  =  1 .N 

BO  =  BO  -  1.0 
C  =  C»(B0/( AO  +  BO) ) 

U  =  ALOG(C)  +  U 

70  Z  =  (A»LNX  +  B*LNY)  -  U 
BO  =  BO  -  1.0 

APB  =  AO  +  BO 
IF  (APB  .GT.  1 .0)  GO  TO  71 
T  =  1.0  +  GAMKAPB  ) 

GO  TO  72 

71  U  =  DBLE(AO)  +  DBLE(BO)  -  1.00 
T  =  (1.0+  GAMKU)  )/APB 

72  BRCOMP  =  AO»EXP( Z )• { 1 .0  +  GAM1(B0))/T 
RETURN 

ALGORITHM  FOR  BO  . GE .  8 

80  U  =  GAMLNKAO)  +  ALGO  I  V(  AO.  BO ) 

Z  =  (A+LNX  +  B'LNY)  -  U 
BRCOMP  =  AO»EXP(Z) 

RETURN 


PROCEDURE  FOR  A  .GE.  8  AND  B  .GE.  8 


100  IF  (A  .GT.  B)  GO  TO  101 

H  *  A/B 

XO  »  H/( 1 .0  +  H) 

YO  =  1 .0/( 1 .0  +  H) 

LAMBDA  =  A  -  (A  +  BI'X 
GO  TO  110 

101  H  »  B/A 

XO  =  1 .0/( 1 .0  +  H) 

YO  =  H/( 1 .0  +  H) 

LAMBDA  *  (A  +  B)>Y  -  B 

1 10  E  =  -LAMBDA/A 

IF  (AB5(E)  .GT.  0.61  GO  TO  HI 
U  =  RLOGKE  ) 

GO  TO  120 

111  U  =  E  -  ALOG(X,/XO) 

120  E  =  LAMBDA/B 

IF  (ABS(E)  .GT.  0.6)  GO  TO  121 
V  =  RL0G1 ( E ) 

GO  TO  130 

121  V  =  E  -  ALOG(Y/YO) 

130  Z  =  EXP( - ( A*U  +  B*V ) ) 

BRCOMP  =  CONST-SORT{B-XO)*Z*EXP( -BCORR(A,B) ) 

RETURN 

END 
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SUBROUTINE  BGRAT(A,  B.  X.  Y.  W.  EPS.  lERR) 


C . . . . — . - . 

C  ASYMPTOTIC  EXPANSION  FOR  IX(A.B)  WHEN  A  IS  LARGER  THAN  B. 

C  THE  RESULT  OF  THE  EXPANSION  IS  ADDED  TO  W.  IT  IS  ASSUMED 

C  THAT  A  .GE.  15  AND  B  . LE .  1.  EPS  IS  THE  TOLERANCE  USED. 

C  lERR  IS  A  VARIABLE  THAT  REPORTS  THE  STATUS  OF  THE  RESULTS 

C . . . . . . 

REAL  J,  L,  LNX.  NU .  N2 
REAL  C(30) .  D(30) 

C 

BM1  =  (B  -  0.5)  -  0.5 
NU  =  A  +  0.5*BM1 
IF  (Y  .GT.  0.375)  GO  TO  10 
LNX  =  ALNREL(-Y) 

GO  TO  11 

10  LNX  =  ALOG(X) 

112=  -NU*LNX 

IF  (B'Z  .EO.  0.0)  GO  TO  100 

computation  of  the  expansion 

SET  R  =  EXP( -Z)»Z*»B/GAMMA(B) 

R  =  B'(1.0  +  GAM1 ( B ) )*EXP(B»AL0G( Z ) ) 

R  =  R»EXP( A»LNX ) ‘EXPIO. 5»BM1»LNX) 

U  =  ALGDIV(B,A)  +  B'ALOG(NU) 

U  »  R*EXP(-U) 

IF  (U  .EO.  0.0)  GO  TO  100 
CALL  GRATKB.Z.R. P.0, EPS) 

C 

V  =  0.25-( t.0/NU)*-2 
T2  =  0.25*LNX*LNX 
L  =  W/U 
U  =  0/R 
SUM  =  d 
T  =  1 .0 
CN  =  1.0 
N2  =  0.0 
DO  22  N  =  1,30 
BP2N  =  e  +  N2 

J  =  (BP2N«(BP2N  +  1.0)*J  +  (Z  +  BP2N  +  1,0)*T)»V 
N2  =  N2  +  2.0 
T  =  T.T2 

CN  =  CN/(N2*lN2  +  1.0)) 

C(N)  =  CN 
S  =  0.0 

IF  (N  .EO.  1 )  GO  TO  21 
NM1  =  N  -  1 
COEF  =  B  -  N 
DO  20  I  =  1 , NM 1 

S  =  S  +  COEF*C( I )*D(N-I ) 

20  COEF  =  COEF  +  B 

21  D(N)  =  BMI'CN  +  S/N 
Dd  =  D(N)*d 

SUM  =  SUM  +  Dd 

IF  (SUM  .LE.  0.0)  GO  TO  100 

IF  (ABS(Dd)  .LE.  EPS*(SUM  +  L))  GO  TO  30 

22  CONTINUE 
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ADD  THE  RESULTS  TO  W 


30  lERR  =  O 

W  =  W  +  U*SUM 
RETURN 


THE  EXPANSION  CANNOT  BE  COMPUTED 


100  I  ERR  =  1 
RETURN 
END 
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SUBROUTINE  GRAT1  ( A , X , R , P , 0 , EPS ) 

REAL  J,  L 

C . - . . . - . 

C  EVALUATION  OF  THE  INCOMPLETE  GAMMA  RATIO  FUNCTIONS 

C  P(A,X)  AND  O(A.X) 

C 

C  IT  IS  ASSUMED  THAT  A  . LE .  1.  EPS  IS  THE  TOLERANCE  TO  BE  USED. 

C  THE  INPUT  ARGUMENT  R  HAS  THE  VALUE  E -X ) ‘X »• A/GAMMA ( A ) . 

C . - . - . - . 

IF  (A»X  .EO.  0.0)  GO  TO  130 
IF  (A  .EO.  0.5)  GO  TO  120 
IF  (X  . LT.  1 . 1)  GO  TO  10 
GO  TO  50 

TAYLOR  SERIES  FOR  P(A.X)/X»*A 

10  AN  =  3.0 
C  =  X 

SUM  =  X/(A  +  3.0) 

TQL  =  0. 1*EPS/( A  +  1 .0) 

11  AN  =  AN  +  1 .0 
C  =  -C.(X/AN) 

T  =  C/(A  +  AN) 

SUM  =  SUM  +  T 

IF  (ABS(T)  .GT.  TOL)  GO  TO  11 
U  *  A«X»( (SUM/6.0  -  0.5/(A  +  2.0))«X  +  1.0/(A  +  1.0)) 

2  =  A»ALOG(X) 

H  =  GAMI(A) 

G  «  1 .0  +  H 

IF  (X  .LT.  0.25)  GO  TO  20 

IF  (A  .LT.  X/2.59)  GO  TO  40 
GO  TO  30 

20  IF  (2  .GT.  -.13394)  GO  TO  40 

30  W  *  EXP( 2  ) 

P  -  WG»(0.5  +  (0.5  -  J)) 

0  =  0.5  +  (0.5  -  P) 

RETURN 

40  L  =  REXP(2) 

W  =  0.5  +  (0.5  +  L) 

0  =  (W*J  -  L)*G  -  H 
IF  (0  . LT .  0.0)  GO  TO  110 
P  =  0.5  +  (0.5  -  0) 

RETURN 

CONTINUED  FRACTION  EXPANSION 

50  A2NM1  =  1.0 
A2N  »  1.0 
B2NM1  =  X 

B2N  =  X  +  ( 1 .0  -  A) 

C  *  1.0 

51  A2NM1  =  X*A2N  +  C*A2NM1 
B2NM1  =  X*e2N  +  C»B2NM1 
AMO  =  A2NM1/B2NM1 
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c  =  c  +  1.0 

CMA  =  C  -  A 
A2N  =  A2NM1  +  CMA*A2N 
B2N  =  B2NM1  +  CMA*B2N 
ANO  =  A2N/B2N 

IF  (ABS(ANO  -  AMO)  . GE .  EPS»ANO)  GO  TO  51 
0  =  R»ANO 

P  =  0.5  +  (0.5  -  0) 

RETURN 

SPECIAL  CASES 

100  P  =  0.0 
0=1.0 
RETURN 

1 10  P  =  1.0 
0  =  0.0 
RETURN 

120  IF  (X  .GE.  0.25)  GO  TO  121 
P  =  ERF(SORT(X)) 

0  =  0.5  +  (0.5  -  P) 

RETURN 

121  0  =  ERFC1 (O.SORT(X) ) 

P  =  0.5  +  (0.5  -  0) 

RETURN 

130  IF  (X  .LE.  A)  GO  TO  100 
GO  TO  1  10 
END 
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REAL  FUNCTION  BASYM(A,  B.  X.  Y.  LAMBDA,  EPS) 


ASYMPTOTIC  EXPANSION  FOR  IX(A.B)  FOR  LARGE  A  AND  B. 
LAMBDA  =  (A  +  B)*Y  -  B  AND  EPS  IS  THE  TOLERANCE  USED. 
IT  IS  ASSUMED  THAT  LAMBDA  IS  NONNEGATIVE  AND  THAT 
A  AND  B  ARE  GREATER  THAN  OR  EQUAL  TO  15. 


REAL  JO.  J1.  LAMBDA 

REAL  AO(21),  BO(21).  C(21).  D(21) 


nUM  IS  THE  MAXIMUM  VALUE  THAT  N  CAN  TAKE  IN  THE  DO  LOOP 
ENDING  AT  STATEMENT  50.  IT  IS  REQUIRED  THAT  NUM  BE  EVEN. 
THE  ARRAYS  AO,  BO.  C,  D  HAVE  DIMENSION  NUM  +  1. 

DATA  NUM/20/ 


EO  =  2/S0RT(PI) 
E1  =  2*»(-3/2) 


DATA  EO/1 . 12837916709551/.  E 1 /. 353553390593274/ 


BASYM  =  0.0 
IF  (A  .GE.  B)  GO  TO  10 
H  =  A/B 

RO  *  1 .0/( 1 .0  +  H) 

R1  ■=  (B  -  A)/B 

WO  '  1 .0/SORT< A-( 1.0  +  H) ) 

GO  TO  20 
10  H  «  B/A 

RO  =  1 . 0/ ( 1 . 0  +  H ) 

R1  »  (B  -  A)/A 

WO  =  1 .0/S0RT(B»( 1 .0  +  H)) 

C 

20  F  =  A»RL0G1 ( -LAMBOA/A )  +  B*RLOG 1 ( LAMBDA/B ) 
T  =  EXP(-F) 

IF  (T  .EO.  0.0)  RETURN 
ZO  =  SQRTIF  ) 

Z  =  0.5»(Z0/E1 ) 

Z2  =  F  +  F 
C 

A0( 1 )  «  (2.0/3.0)»R1 
C(  1 )  =  -  0. 5*A0(  1  ) 

0(1)  =  -  C(  1  ) 

JO  =  (0.5/E0)*ERFC1( 1 ,Z0) 

J1  =  El 

SUM  *  JO  +  D( 1 )«W0*J1 
C 

S  »  1 .0 
H2  *  H*H 
HN  =  1.0 
W  =  WO 
ZNM1  =  Z 
ZN  «  Z2 

DO  50  N  =  2,  NUM.  2 
HN  =  H2»HN 

AO(N)  =  2.0*R0*(1.0  +  H«HN)/(N  +  2.0) 
NP1  =  N  +  1 
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S  =  S  +  HN 

AO(NPI)  =  2.0*R1«S/(N  +  3.0) 

00  41  I  »  W,  UP^ 

R  =  -0.5-( 1  +  1.0) 

B0(  1 )  -  R*A0(  1  ) 

00  31  M  «  2.  I 
BSUM  «  0.0 
MM1  »  M  -  1 
00  30  J  •  1 .  MM 1 
MMo  -  M  -  u 

30  BSUM  =  BSUM  ♦  (J«R  -  MMd ) -A0( J ) •B0( MMU ) 

31  BO(M)  •  R*A0(M)  +  eSUM/M 
C( I )  -  B0( I  )/( 1  ♦  1.0) 

DSUM  •  0.0 
IM1  «  I  -  1 

00  40  J  =  1  .  I M 1 
IMJ  »  I  -  d 

40  OSUM  =  DSUM  +  D(IMd)*C(o) 

41  D( I  )  »  -(OSUM  +  C( I  )  ) 

dO  •  E1»ZNM1  +  (N  -  1.0)»d0 

d1  •  E 1 'ZN  +  N*d1 

ZNM1  «  22-ZNM1 

ZN  «  Z2*ZN 

W  *  WO"W 

TO  «  D(N)*W*dO 

W  «  WO  *  w 

T1  •  0(NP1  )«W*d1 
SUM  «  SUM  +  (TO  +  T1 ) 

IF  ((ABS(TO)  +  ABS(TI))  .LE.  EPS'SUM)  go  TO  60 
50  CONTINUE 

60  U  »  EXP( -BCORR( A.B)  ) 

BASYM  *  E0*T*U*SUM 

RETURN 

END 
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REAL  FUNCTION  REXP(X) 


COMPUTATION  OF  EXP(X)  -  1 


DATA  P1/  .914041914819518E-09/.  P2/  . 23S082361044469E-01 / . 
01/-.4999999990859S8E+00/.  02/  .  107  14  156a980644E-t-O0/ , 
03/- . 1 19041 179760821E-01/.  04/  . 595 1 308 1 1860248E -03/ 


IF  (ABS(X)  .GT.  0.15)  GO  TO  10 

REXP  «  X»(((P2»X  +  P1)*X  +  1 .0)/( ( ( (04»X  ♦  03)*X  +  02 ) *X 
•  ♦  01 )»X  +  1.0)) 

RETURN 

C 

10  W  -  EXP(X) 

IF  (X  .GT.  0.0)  GO  TO  20 
REXP  •  (W  -  0.5)  -  0.5 
RETURN 

20  REXP  «  W»(O.S  +  (0.5  -  1.0/W)) 

RETURN 

END 
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REAL  FUNCTION  ALNREL(A) 


EVALUATION  OF  THE  FUNCTION  LN( 1  +  A) 


DATA  P1/- . 129418923021993E+01/.  P2/ . 405303492862024E+00/ . 

•  P3/- . 178874546012214E-01/ 

DATA  01/- . 162752256355323E+01/.  02/ . 7478 1 10140376 16E+00/ , 

•  03/- .845104217945565E-01/ 

C . 

IF  (ABS(A)  .GT.  0.375)  GO  TO  10 
T  =  A/(A  +  2.0) 

T2  =  T*T 

W  =  (((P3»T2  +  P2)»T2  +  P1)«T2  +  1.0)/ 

•  ({(03*T2  +  02)*T2  +  01)»T2  +  1.0) 

ALNREL  =  2.0»T*W 

RETURN 

C 

10  X  =  1 .DO  +  DBLE(A) 

ALNREL  =  ALOG(X) 

RETURN 

END 
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REAL  FUNCTION  RLOGKX) 


COMPUTATION  OF  X  -  LN{1  +  X) 


DATA  A/.566749439387324E-01/ 
DATA  B/.4565126O8815524E-01/ 


DATA  PO/  . 333333333333333E+00/ .  P 1 /-. 2246964 1 3 1 1 2536E+00/ . 
•  P2/  .620886815375787E-02/ 

DATA  01/- . 1274089239336238+01/.  02/  . 3545087 1 8369557E+00/ 


IF  (X  .LT.  -0.39  .OR.  X  . GT .  0.57)  GO  TO  100 
IF  (X  . LT.  -0.18)  GO  TO  10 
IF  (X  .GT.  0. 18)  GO  TO  20 

ARGUMENT  REDUCTION 


H  =  X 
W1  =  0.0 
GO  TO  30 

10  H  =  DBLE(X)  +  0.3D0 
H  =  H/0.7 
W1  =  A  -  H*0.3 
GO  TO  30 

20  H  =  0.75D0*DBLE(X)  -  0.25D0 
W1  =  B  +  H/3.0 

SERIES  EXPANSION 

30  R  »  H/(H  +  2.0) 

T  =  R«R 

W  =  ((P2»T  +  P1)«T  +  P0)/((02*T  +  01)»T  +  1.0) 
RL0G1  =  2.0»T«( 1.0/( 1.0  -  R)  -  R*W)  +  W1 
RETURN 


100  W  =  (X  +  0.5)  +  0.5 
RL0G1  =  X  -  ALOG(W) 
RETURN 
END 
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REAL  FUNCTION  ERF  (X) 

EVALUATION  OF  THE  REAL  ERROR  FUNCTION 


C 


C 


C 


C 


C 


REAL  A(4),B(4).P(8).0(8).R(5).S(5) 
DATA  C/. 564189583547756/ 


DATA 

DATA 

DATA 

DATA 

DATA 

DATA 


A( 1 )/-1 .65581836870402E-4/ 
A( 3 )/ 1 .02201 1369 18406E-1/, 
B{ 1 )/4.64988945913179E-3/. 
B(3)/4.239067326832O1E-1/. 
P( 1 )/-1 .36a64857382717E-7/ 
P(3)/7.21 175825088309E00/, 
P(5)/1 .52989285046940E02/. 
P( 7 )/4. 5 191895371 1873802/, 
0( 1 ) / 1 . OOOOOOOOOOOOOOEOO/ . 
0( 3 )/7 . 70001529352295E01/, 
0(5 )/6. 3898026446563 1E02/. 
0( 7 ) /7 . 90950925327e98E02/ . 
R( 1  )/2.  10144 126479064E00/ , 
R(3)/2. 13688200555087E01/ . 
R( 5 )/2. 8209479 1773523E-1/ 
S( 1 )/9.41537750555460E01/. 
S( 3 )/9. 9019 18 146239 14E01/. 
S ( 5  )  / 1 . OOOOOOOOOOOOOOEOO/ 


A( 2 )/3. 253240983577386 -2/. 
A(4)/1  .  12837916709552E00/ 
B( 2 )/7. 013334 17 15851  IE -2/. 
B ( 4 ) / 1 . OOOOOOOOOOOOOOEOO/ 
P(2)/5.64195517478974E-1/. 
P(4)/4 . 31622272220567E01/. 
P( 6 ) /3 . 393208 16734344E02/ . 
P( 8 )/3. 0045926 10201 62E02/ 
0(2>/1 . 27827273 196294E01/. 
0( 4  )/2 . 77585444  743988602/ . 
0( 6 )/9 . 3 1 3540948506 10E02/ . 
0( 8 )/3 . 00459260956983E02/ 
R( 2 )/2. 6237014 1675 169E01/, 
R( 4 )/4 . 658078287 18470EOO/ , 

S(2)/1 .87114811799590E02/, 
S(4)/1 .80124575948747E01/. 


AX  =  ABS(X) 

IF  (AX  .GE.  0.5)  GO  TO  10 
T  =  X*X 

TOP  =  ((A(1).T  +  A(2))*T  +  A(3))*T  +  A(4) 
80T  =  ((8(1)*T  +  B(2))*T  ♦  B(3))*T  +  B(4) 
ERF  «  X*TOP/BOT 
RETURN 

10  IF  (AX  .GT.  4.0)  GO  TO  20 


TOP 

a 

(((((( P( 1 ) -AX 

-► 

P(2) )-AX 

+ 

P(  3) 

)-AX 

P(4 ) )-AX 

+  P(5) )-AX 

• 

+ 

P<6 ) )-AX 

+ 

P(7) 

)-AX 

+ 

P(8) 

BOT 

a 

( ( ( ( ( (0( 1 )-AX 

+ 

0( 2 ) )-AX 

+ 

0(3) 

)  -AX 

+ 

0( 4 ) ) -AX 

+  0(5)) -AX 

0 

+ 

0(6) )-AX 

+ 

0(7) 

)»AX 

•f 

0(8) 

ERF 

= 

0.5  +  (0.5  - 

EXP(  -X-X)-TOP/BOT) 

IF 

(X 

.LT.  0.0)  ERF 

= 

-ERF 

RETURN 

20  ERF  =  1.0 

IF  (AX  .GE.  5.6)  GO  TO  21 
X2  =  X*X 
T  =  1.0/X2 

TOP  =  (((R(1)*T  +  R(2))-T  *  R(3))»T  +  R(4))»T  +  R(5) 
BOT  =  (((S(1)*T  +  S(2))*T  +  S(3))-T  +  S(4))»T  +  S(5) 
ERF  =  (C  -  TOP/( X2»B0T)  )  /  AX 
ERF  =  0.5  +  (0.5  -  EXP(-X2)*ERF) 

21  IF  (X  .LT.  0.0)  ERF  =  -ERF 
RETURN 

END 
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REAL  FUNCTION  ERFC1  (IND,  X) 

EVALUATION  OF  THE  COMPLEMENTARY  ERROR  FUNCTION 

ERFCKIND.X)  =  ERFC(X)  IF  IND  =  0 

ERFCKIND.X)  =  EXP(X»X)»ERFC(X)  OTHERWISE 


REAL  A(4),B(4),P(8),0(8),R(5).S<5) 
DOUBLE  PRECISION  W 


DATA  C/. 564189583547756/ 


DATA  A( 1 )/- 1 .6S581836870402E-4/. 

A( 3 )/ 1 .02201 13691 8406E- 1/ . 
DATA  B( 1 )/4.64988945913179E-3/. 

B( 3 )/4 . 23906732683201 E- 1/ . 
DATA  P( 1 )/-1 .36864857382717E-7/, 
P(3)/7.21 17582508e309E00/. 
P(5)/1 .52989285046940E02/. 
P(7)/4.51918953711873E02/, 
DATA  0( 1 )/ 1  - OOOOOOOOOOOOOOEOO/ . 
0( 3)/7 . 70001529352295E01/, 
0(5 )/6. 3898026446563 1E02/. 
0(7 )/7 .90950925327898E02/, 
DATA  R(  1  )/2 .  10144 126479064EOO/ . 
R(3)/2. 13688200555087E01/. 
R( 5 )/2. 8209479 1773523E-1/ 
DATA  S( 1 )/9.41537750555460E01/. 
S( 3 )/9. 9019 181 46239 14E01/. 

S ( 5 ) / 1 . OOOOOOOOOOOOOOEOO/ 


A(2)/3. 25324098357738E-2/ , 
A(4)/1 . 128379 16709552E00/ 
B( 2 )/7. 013334 17 15851  IE -2/. 
B( 4 )/ 1 . OOOOOOOOOOOOOOEOO/ 
P(2)/5.64195517478974E- 1/, 
P(4  )/4. 31622272220567E01/, 
P( 6 )/3. 393208 16734344E02/, 
P( 8 )/3. 0045926 1020162E02/ 
0(2)/1 . 27827273 196294 E01/, 
0( 4 )/2 . 77585444743988E02/ , 
0( 6 )/9. 31 3540948506 10E02/. 
0( 8 )/3 . 00459260956983E02/ 
R( 2 )/2 .6237014 1675 169E01/ . 
R( 4 )/4 . 658078287 1 8470E00/ , 

S( 2)/1 .871 1481 1799590E02/. 
S(4)/1 .80124575948747E01/. 


ABS(X)  .LT.  0.47 


AX  ==  ABS(X) 

IF  (AX  .GE.  0.47 )  GO  TO  10 
T  =  X»X 

TOP  =((A(1)«T  +  A(2))*T  +  A(3))-T  +  A(4) 
BOT  =  ((B(1)'T  +  B(2))-T  +  B(3))*T  +  B( 4 ) 
ERFC1  =  0.5  +  (0.5  -  X*TOP/BOT) 

IF  (IND  .NE.  0)  ERFC1  =  EXP(T)  •  ERFC1 
RETURN 


0.47  .LE.  ABS(X)  .LE.  4 


C 


10  IF  (AX 

.GT.  4.0)  GO 

TC 

1  30 

TOP  = 

(((((( P(  1  )*AX 

+ 

P(  2 

» 

•¥ 

P(6 

BOT  = 

( ( ( ( ( (0( 1 )»AX 

4 

0(2 

ERFC1 

*  TOP/BOT 

4 

0(6 

)  )»AX 

+ 

P( 3) )«AX 

4 

P(4  ) 

)*AX 

)  )*AX 

4 

P(7) )*AX 

4 

P(8) 

)  J-AX 

4 

0(3)  )»AX 

4 

0(4) 

)»AX 

)  )»AX 

4 

0(7) ).AX 

4 

0(8) 

+  P(5))»AX 
+  0(5))*AX 


20  IF  (IND  .EO.  0.0)  GO  TO  21 

IF  (X  .LT.  0.0)  ERFC1  =  2.0*EXP(X*X)  -  ERFC1 
RETURN 

2 1  W  =  DBLE ( X )*DBLE ( X  ) 

T  =  W 
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E  =  W  -  DBLE(T) 

ERFC1  =  ((0.5  +  (0.5  -  E>)  •  EXP(-T))  »  ERFC1 

IF  (X  .LT.  0.0)  ERFC1  =  2.0  -  ERFC1 

RETURN 


ABS(X)  ,GT.  4 

30  IF  (X  .LE.  -5.5)  GO  TO  40 

IF  (IND  .EO.  0.0  .AND.  X  .GT.  50.0)  GO  TO  50 
T  =  ( 1 ,0/X)*»2 

TOP  =  (((R(1)*T  +  R(2))»T  +  R(3))»T  +  R(4))»T  +  R(5) 
BOT  =  (((S(1)»T  +  S(2))«T  +  S(3))»T  +  S(4))*T  +  S(5) 
ERFC1  =  (C  -  T»TOP/BOT)  /  AX 
GO  TO  20 


LIMIT  VALUE  FOR  LARGE  NEGATIVE  X 
40  ERFC1  =  2.0 

IF  (INO  .NE.  0)  ERFC1  =  2.0*EXP(X«X) 
RETURN 

LIMIT  VALUE  FOR  LARGE  POSITIVE  X 
WHEN  IND  =  O 

50  ERFC1  =  0.0 
RETURN 
END 
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real  function  GAMI(A) 


computation  of  1/GAMMA(A+1)  -  1  FOR  -0.5  . LE .  A  .LE.  1.5 


C 

C 

C 

C 

c 

c 

c 


REAL 

P(  7  )  , 

0( 5  )  .  R( 9) 

DATA 

)* 

P(  1  )  / 
P(3)/- 
P  (  5  ) , ' 
P(  7  )/ 

. 5772 1 566490 1533E+00/. 

. 2309753808576755+00/ , 

. 766968 18 1649490E -02/ , 
.58959742861 1429E-03/ 

P( 2 )/- 

P(4)/ 

P(6)/- 

.  409078 193005776E+00/ 

.  597275330452234E-01/ 
.5  1488977 1323592 E -02/ 

DATA 

X. 

0(1)7 

0(3)/ 

0(5)/ 

.  1 0OOOOOOOOOOOOOE  +01/, 

. 158451672430138E+00/, 
.4  23244 29789696  IE -02/ 

0(2)/ 

0(4)/ 

.  4275696130952 14E  +  00/ 
.261 132021441447E-01/ 

DATA 

R(  1  )/- 
R( 3  ;/- 
R(5  1/ 

R  (  7  \  / 
R( 9  )/- 

.  4227843350984688+00/ , 

. 244757765222226E+00/ , 

. 930357293360349E-03/. 
.223047661 158249E-02/, 

. 132674909766242E-03/ 

R(2 )/- 
R(4  )/ 
R(6)/- 
R(8)/ 

.  771330383816272E+00/ 

.  1  18378989872749E+00/ 

.  1  18290993445146E-01/ 

.  266505979058923E -03/ 

DATA 

S1  / 

. 273076 135303957E+00/. 

S2  / 

.  559398236957378E-01/ 

T  *  A 

0  =  A  -  0.5 

IF  (D  .O'"  0.0)  T  =  D  -  0.5 

IF  (7 )  30, 10,20 

10  GAMl  =  0.0 
RETURN 

20  TOP  =  (l(((P(7)»T  +  P(6))»T  +  P(5))*T  +  P(4))-T  +  P ( 3  )  )  • T 

•  +  P(  2  )  )*T  +  P(  1  ) 

EOT  =  (((0(5)*T  *  0(4))-T  +  0(3))'T  +  0(2))*T  +  1.0 
W  =  T0P/60T 

IF  (D  .GT.  0.0)  GO  TO  21 
GAMl  -  A*W 
RETURN 

21  GAMl  -  (T/A)'((w  -  0.5)  -  0.5) 

RETURN 

30  TOP  =  ((((((( R( 9  ) -T  +  R(8))*T  +  R(7))*T  *  R(6))-T  +  R(5))‘T 

+  R(4))*T  +  R(3))*T  +  R(2))'T  +  R(1) 

EOT  =  ( S2'T  +  S 1  )  -T  +  1.0 
W  =  TOP/eOT 

IF  (0  .GT.  0.0)  GO  TO  31 

GAM  1  =  A*((W  0.5)  +  0.5) 

RETURN 

31  GAMl  =  T*W/4 
RETURN 

END 
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REAL  FUNCTION  GAMLN(A) 

C  evaluation  of  LN(GAMMA(A))  for  positive  a 

C  WRITTEN  BY  ALFRED  H.  MORRIS 

C  NAVAL  SURFACE  WEAPONS  CENTER 

C  DAHLGREN.  VIRGINIA 

C  - - - 

C  D  =  0.5-( LN(2»PI )  -  1 ) 

C  - 

DATA  0/ . 4 18938533204673/ 

C  - 

DATA  00/ . 833333333333333E-01/ .  C 1 /-. 27777777777048 1 E -02/ , 

1  C2/ . 793650663183693E-03/ .  C3/ 595 1 563364 2859 1 E -03/ . 

2  C4/  820756370353826E-03/ 

C  ---  - - - - - - - 

IF  (A  . GT .  0.8)  GO  TO  10 

GAMLN  -  GAMLNKA)  -  ALOG(  A  ) 

RETURN 

10  IF  I  A  .GT .  2. 25  )  GO  TO  20 
X  =  DBLE(A)  -  1.D0 

GAMLN  =  GAMLNKX) 

RE  TURN 

r 

20  IF  (A  .GE.  15.0)  GO  TO  30 

N  =  A  •  1  25 

X  =  A 
W  =  1.0 
DO  21  I  '  1.N 
X  =  X  •  1.0 

21  W  =  X-W 

GAMLN  =  GAMLNKX  -  1.0)  +  ALOG(W) 

RETURN 

C 

30  X  =  (  1  .O^A  )-'2 

W  -  l((:C4»X  +  C3)-X  +  C2)*X  +  C1)*X  +  C0)/A 
GAMLN  =  (D  +  W)  (A  -  0 . 5 ) '  1  ALOG(  A  )  -  1.0) 

END 
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REAL  FUNCTION  GAMLNKA) 


EVALUATION  OF  LN(GAMMA(1  +  A))  FOR  -0.2  .LE.  A  . LE .  1.25 


DATA  PO/  . 577215664901533E+00/.  Pi/  . 844203922 187225E+00/ . 
P2/- . 168860593646662E+O0/.  P3/- . 7804276 1553359 1 E+OO/ . 
P4/-. 4020557993 10489E+00/.  P5/ -. 6735622 1 432567 1 E -0 1 / , 
P6/- .271935708322958E-O2/ 

DATA  01/  . 288743195473681E+01/.  02/  . 3 1 27550889 1 4843E +0 1 / . 
03/  .1568751932950396+01/.  04/  . 36 1 95 1 990 101 499E+00/ . 
05/  . 325038868253937E-01/ .  06/  . 6674656 1 8796 1 64E -03/ 


DATA  RO/ . 422784335098467E+00/.  R 1 /. 8480446 1 4534529E+00/ . 
R2/ .56522 105069 1933E+00/.  R3/ . 1 565 1 306048655 1 E +0O/ . 

R4/. 1705024840226506-01/.  R5/ . 497958207639485E -03/ 

DATA  SI/. 124313399877507E+O1/.  52/ . 548042 109832463E+0O/ , 

53/ . 101552187439830E+00/ .  54/ . 7 1 33096 1 239 1 OOOE -02/ . 

55/ . 1 161654759896 16E-03/ 


IF  (A  .GE.  0.6)  GO  TO  10 

W  =  ((((((P6»A  +  P5)»A  +  P4)*A  +  P3)*A  +  P2>»A  +  P1)»A  +  PO)/ 
■  ((((((06*A  +  05)»A  +  04)-A  +  03)»A  +  02)»A  +  01)'A  +  1.0) 

GAMLN1  =  -A»W 
RETURN 
C 

10  X  *  DBLE( A)  -  1 .DO 

W  =  (((((R5»X  +  R4)*X  +  R3)»X  ♦  R2)*X  +  R1)*X  +  RO)/ 

•  (((((S5*X  +  S4)*X  +  S3)»X  +  S2)*X  +  S1)*X  +  1.0) 

GAMLN1  =  X»W 

RETURN 

END 
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REAL  FUNCTION  BETALN(AO.  BO) 


EVALUATION  OF  THE  LOGARITHM  OF  THE  BETA  FUNCTION 


E  =  0.5»LN(2»PI) 

DATA  E  /. 918938533204673/ 


A  =  AMINKAO.BO) 

B  =  AMAXKAO.  :0) 

IF  (A  .GE.  8.0)  GO  TO  60 
IF  (A  .GE.  1.0)  GO  TO  20 


PROCEDURE  WHEN  A  . LT .  1 


IF  (B  .GE.  8.0)  GO  TO  10 

BETALN  =  GAMLN(A)  +  (GAMLN(B)  -  GAMLN(A  +  B)) 
RETURN 

10  BETALN  =  GAMLN(A)  +  ALGDIV(A.B) 

RETURN 


PROCEDURE  WHEN  1  . LE .  A  .  LT  .  8 


20  IF  (A  .GT.  2.0)  GO  TO  30 
IF  (B  .GT.  2.0)  GO  TO  21 

BETALN  =  GAMLN(A)  +  GAMLN(B)  -  G5UMLN(A,B) 
RETURN 


IF  (8  . LT.  a. 0)  GO  TO  40 

BETALN  =  GAMLN(A)  +  ALGDIV(A,B) 
RETURN 


REDUCTION  OF  A  WHEN  A  .GT.  2 


30  N  =  A  -  1.0 
W  =  0.0 

DO  31  I  =  1  .N 
A  =  A  -  1 . 0 

H  =  A/B 

W  =  W  +  ALOG(H/( 1.0  +  H) ) 

31  CONTINUE 

IF  (B  LT.  8,0)  GO  TO  40 

BETALN  =  W  +  GAMLN(A)  +  ALGDIV(A,B) 

RETURN 


REDUCTION  OF  B  WHEN  B  .LT.  8 


40  N  =  B  -  1.0 
Z  =  1.0 

DO  4  1  I  =  1  .N 
B  =  B  -  1 .0 

41  Z  =  Z*(8/(A  +  B)) 

BETALN  =  W  +  ALOG(Z)  +  (GAMLN(A)  +  (GAMLN(B)  -  GSUMLN( A , B ) ) ) 
RETURN 


PROCEDURE  WHEN  A  . GE .  8 


0-29 


60  W  =  BC0RR(A,B) 

H  =  A/B 

C  =  H/  (  1 . 0  +  H  ) 

U  =  -(A  -  O.5)*AL0G(C) 

V  =  B»AUNREL(H) 

IF  (U  . LE.  V)  GO  TO  61 

BETALN  =  ( ( ( -0.5*AL0G(B)  +  E)  +  W)  - 
RETURN 

61  BETALN  =  ( ( ( -0. 5»AL0G( B )  +  E)  +  W)  -  U) 
RETURN 

END 


V)  -  U 
-  V 
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REAL  FUNCTION  GSUMLN(A.B) 


EVALUATION  OF  THE  FUNCTION  LN(GAMMA(A  +  B)) 
FOR  1  .LE.  A  .LE.  2  AND  1  . LE .  B  . LE .  2 


X  =  DBLE(A)  +  DBLE(B)  -  2. DO 
IF  (X  .GT.  0.25)  GO  TO  10 
GSUMLN  =■  GAMLN  1(1.0  +  X) 

RETURN 

10  IF  (X  .GT.  1.25)  GO  TO  20 

GSUMLN  =  GAMLNKX)  +  ALNREL(X) 

RETURN 

20  GSUMLN  «  GAMLNKX  -  1.0)  +  AL0G(X»(1.0  +  X)) 
RETURN 
END 
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REAL  FUNCTION  BCORR(AO.  SO) 


EVALUATION  OF  DEL(AO)  +  DEL(BO)  -  DEL(AO  <■  BO)  WHERE 
LNrQAMMA(A))  »  (A  -  0.5)-LN(A)  -  A  +  C.5-LN(2*PI)  +  OEL(A). 
IT  IS  ASSUMED  THAT  AO  . GE .  8  AND  BO  . GF ,  8. 


DATA  CO/.833333333333333E-01/.  Cl . 2‘'777777776099 'E-02/ . 
C2/.793650666825390E-03/.  C3/-  59520293 1 35 1870E -03/ , 
C4/ .83730803403 12 15E -03/ .  C5/- . 1653229627807 1 3E -02/ 


A  =  AMINKAO,  BO) 

B  -  AMAXKAO.  BO) 

H  =  A/B 

C  =  H/ (  1 . 0  +  H ) 

X  =  1.0/ (1.0+  H) 

X2  »  X*X 

SET  SN  *  ( 1  -  X»-N)/{ 1  -  X ) 


S3  «  1.0+ 

(X 

> 

X2  ) 

S5  •  1 .0  + 

(X 

X2»S3  ) 

S7  =  1 .0  + 

(X 

X2-S5) 

S9  *  1 .0  + 

(  X 

X2-S7  ) 

SI  1  =  1.0  + 

(X 

+  X2*S9) 

SET  W  *  DEL(B)  -  DELIA  +  B) 

T  »  (1.0/B)«*2 

W  *  ((((CS»S11*T  +  C4«S9)*T  +  C3»S7)-T  +  C2-S5)*T  +  C1*S3)-T  +  CO 
W  «  W(C/B) 


COMPUTE  DELIA)  +  W 


T  =  ( 1 .0/A)*-2 

BCORR  =  (((((C5*T  +  C4)-T  +  C3)«T  +  C2)-T  +  C1)»T  +  C0)/A  +  W 

RETURN 

END 
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REAL  FUNCTION  ALGDIV(A.B) 


C 

C 

C 

C 

C 

C 


C 


COMPUTATION  OF  LN( GAMMA { B ) /GAMMA( A+B ) )  WHEN  B  .GE.  8 
IN  THIS  CODE  DEL(Z)  IS  THE  FUNCTION  WHERE 

LN(GAMMA(Z))  =  (Z  -  0.5)»LN(Z)  -  Z  +  0.5»LN(2*P1)  +  DEL(Z) 


DATA  C0/.833333333333333E-01/,  C 1 /-. 27777777776099 1 E -02/ . 

*  C2/.793650666825390E-03/,  C3/- . 59520293 1 35 1870E -03/ , 

*  C4/ .83730803403 12 15E -03/.  C5/- . 1653229627807 1 3E -02/ 


10 


IF 

(A  .LE. 

B)  GO  TO 

10 

H  =  B/A 

C  =  1.0/ 

(1.0 

+  H 

) 

X  =  H/(  1 

.0  + 

H) 

0  =  A  ♦ 

(B  - 

0.5 

) 

GO  TO  20 

H  = 

A/B 

C  = 

H/(  1 .0 

+  H) 

X  = 

1 .0/(1  . 

0  + 

H) 

0  = 

B  +  (A 

-  0. 

5) 

SET 

SN  = 

(  1  - 

X2 

=  X*X 

S3 

=  1.0  + 

(X  + 

X2  ) 

S5 

=  1.0  + 

(X  + 

X2» 

S3) 

S7 

*  1.0  + 

(X  + 

X2» 

S5) 

S9 

=  1.0  + 

(X  + 

X2* 

S7) 

S1 1 

+ 

o 

H 

(X 

+  X2 

•S9) 

SET 

W  = 

DEL(B) 

T  = 

( 1 .0/B) 

••2 

W  = 

( ( ( (C5» 

S11* 

T  + 

C4*S9) 

w  = 

W»(C/B) 

X) 


COMBINE  THE  RESULTS 


U  =  0»ALNREL(A/B ) 

V  =  A*( ALOG(B)  -  1.0) 

IF  (U  .LE.  V)  GO  TO  30 
ALGDIV  =  (W  -  V)  -  U 
RETURN 

30  ALGDIV  =  (W  -  U)  -  V 
RETURN 
END 


0-33 


REAL  FUNCTION  SPMPAR(I) 
INTEGER  I 


C- . - . - . . . - . 

c 

C  SPMPAR  PROVIDES  THE  SINGLE  PRECISION  MACHINE  PARAMETERS  FOR 

C  THE  COMPUTER  BEING  USED.  IT  IS  ASSUMED  THAT  THE  ARGUMENT 

C  I  IS  AN  INTEGER  HAVING  ONE  OF  THE  VALUES  1.  2,  OR  3.  IF  THE 

C  SINGLE  PRECISION  ARITHMETIC  BEING  USED  HAS  T  BASE  B  DIGITS  AND 

C  ITS  SMALLEST  AND  LARGEST  EXPONENTS  ARE  EMIN  AND  EMAX,  THEN 

C 

C  SPMPAR(1)  =  B»»(1  -  T).  THE  MACHINE  PRECISION. 

C 

C  SPMPAR(2)  =  B**{EMIN  -  1).  THE  SMALLEST  MAGNITUDE, 

C 

C  SPMPARO)  =  B»*EMAX*(1  -  B*»(-T)).  THE  LARGEST  MAGNITUDE. 

C 

C  TO  DEFINE  THIS  FUNCTION  FOR  THE  COMPUTER  BEING  USED.  ACTIVATE 

C  THE  data  STATMENTS  FOR  THE  COMPUTER  BY  REMOVING  THE  C  FROM 

C  COLUMN  1.  (ALL  OTHER  DATA  STATEMENTS  IN  SPMPAR  SHOULD  HAVE  C 

C  IN  COLUMN  1 , ) 

C 

C . - . . 

C 

C  SPMPAR  IS  AN  ADAPTATION  OF  THE  FUNCTION  R1MACH.  WRITTEN  BY  P.A. 

C  FOX,  A.D.  HALL.  AND  N.L.  SCHRYER  (BELL  LABORATORIES).  SPMPAR 

C  WAS  DESIGNED  BY  B.S.  GARBOW.  K.E.  HILLSTROM,  AND  J.J.  MORE 

C  (ARGONNE  NATIONAL  LABORATORY).  THE  MAJORITY  OF  PARAMETER  VALUES 

C  ARE  FROM  BELL  LABORATORIES. 

C 

C . - . 

INTEGER  MCHEPS(2) 

INTEGER  MINMAG(2 ) 

INTEGER  MAXMAG(2 ) 

REAL  RMACH(3) 

EQUIVALENCE  ( RMACH( 1 ) . MCHEPS( 1 ) ) 

EQUIVALENCE  ( RMACH( 2 ) . MINMAG( 1 ) ) 

EQUIVALENCE  ( RMACH ( 3 ) . MAXMAG( 1 ) ) 

C 

C  MACHINE  CONSTANTS  FOR  THE  BURROUGHS  1700  SYSTEM. 

C 

C  DATA  RMACH( 1 )  /  Z4EA800000  / 

C  DATA  RMACH(2)  /  Z400800000  / 

C  DATA  RMACH(3)  /  Z5FFFFFFFF  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  BURROUGHS  5700/6700/7700  SYSTEMS. 

C 

C  DATA  RMACH(1)  /  01301000000000000  / 

C  DATA  RMACH(2)  /  01771000000000000  / 

C  DATA  RM4CH(3)  /  00777777777777777  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  CDC  6000/7000  SERIES. 

C  (OCTAL  FORMAT  FOR  FORTRAN  4  COMPILERS) 

C 

C  DATA  RMACH( 1 )  /  1 64 1 4000000000000000B  / 

C  DATA  RMACH(2)  /  000 1 4000000000000000B  / 

C  DATA  RMACH(3)  /  37767777777777777777B  / 

C 
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C  MACHINE  CONSTANTS  FOR  THE  CDC  6000/7000  SERIES. 

C  (INTEGER  FORMAT  FOR  FORTRAN  4  AND  5  COMPILERS) 

C 

DATA  MCHEPS(1)  /  261630990852554752  / 

DATA  MINMAG( 1  )  /  422212465065984  / 

DATA  MAXMAGI 1 )  /  576179277326712831  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  CRAY-1. 

C 

C  data  RMACH( 1 )  /  0377224000000000000000B  / 

C  DATA  RMACH(2)  /  0200034000000000000000B  / 

C  DATA  RMACH(3)  /  0577777777777777777776B  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  DATA  GENERAL  ECLIPSE  S/200. 

C 

C  NOTE  -  IT  MAY  BE  APPROPRIATE  TO  INCLUDE  THE  FOLLOWING  CARD 

C  STATIC  RMACH(3) 

C 

C  DATA  MINMAG/20K,0/.MAXMAG/77777K, 177777K/ 

C  data  MCHEPS/36020K,0/ 

C 

C  MACHINE  CONSTANTS  FOR  THE  HARRIS  220. 

C 

C  DATA  MCHEPS(I)  /  '20000000.  '00000353  / 

C  data  MINMAG(I)  /  '20000000.  '00000201  / 

C  data  MAXMAGC 1 )  /  '37777777.  '00000177  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  HONEYWELL  600/6000  SERIES. 

C 

C  DATA  RMACH( 1 )  /  07 16400000000  / 

C  DATA  RMACH(2)  /  0402400000000  / 

C  DATA  RMACHO)  /  0376777777777  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  HP  2100 

C 

C  DATA  MCHEPS(I).  MCHEPS(2)  /  40000B .  327B  / 

C  DATA  MINMAG(I).  MINMAG(2)  /  40000B .  1  / 

C  DATA  MAXMAGI 1),  MAXMAGI 2 )  /  77777B,  177776B  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  HP  9000 

C 

C  DATA  RMACH( 1 )  /  .1192093E-06  / 

C  DATA  RMACH(2)  /  .5877472E-38  / 

C  DATA  RMACHO)  /  .  3402823E  +  39  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  IBM  360/370  SERIES. 

C  THE  AMDAHL  470/V6.  THE  ICL  2900.  THE  ITEL  AS/6, 

C  THE  XEROX  SIGMA  5/7/9  AND  THE  SEL  SYSTEMS  85/86. 

C 

C  DATA  RMACH( 1)  /  Z3C100000  / 

C  DATA  RMACH(2)  /  Z00100000  / 

C  DATA  RMACHO)  /  Z7FFFFFFF  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  IBM  PC  -  MICROSOFT  FORTRAN 

C 

C  DATA  MCHEPS(I)  /  «34000000  / 

C  DATA  MINMAG( 1 )  /  #00800000  / 
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C  DATA  MAXMAG(I)  /  #7F7FFFFF  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  IBM  PC  -  PROFESSIONAL  FORTRAN. 

C  LAHEY  FORTRAN.  AND  RM  FORTRAN 

C 

C  DATA  MCHEPS(I)  /  Z'34000000'  / 

C  DATA  MINMAG(I)  /  Z' 00800000'  / 

C  DATA  MAXMAG( 1 )  /  Z'7F7FFFFF'  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  PDP-10  (KA  OR  KI  PROCESSOR). 

C 

C  DATA  RMACH( 1)  /  "147400000000  / 

C  DATA  RMACH(2)  /  "000400000000  / 

C  DATA  RMACHO)  /  "377777777777  / 

C 

C  MACHINE  CONSTANTS  FOR  THE  PDP-11  FORTRAN  SUPPORTING 

C  32-BIT  INTEGERS  (EXPRESSED  IN  INTEGER  AND  OCTAL). 

C 

C  DATA  MCHEPS(1)  /  889192448  / 

C  DATA  MINMAG(  1  )  /  8388608  / 

C  DATA  MAXMAG(I)  /  2147483647  / 

C 

C  DATA  RMACH( 1 )  /  006500000000  / 

C  DATA  RMACH(2)  /  000040000000  / 

C  DATA  RMACHO)  /  017777777777  / 

C 

C  MACHINE  constants  FOR  THE  PDP-11  FORTRAN  SUPPORTING 

C  16-BIT  INTEGERS  (EXPRESSED  IN  INTEGER  AND  OCTAL). 

C 

C  DATA  MCHEPS( 1 ).MCHePS(2)  /  13568.  0  / 

C  data  MINMAG( 1 ) .MINMAG(2)  /  128.  0  / 

C  DATA  MAXMAG(  1  ).MAXMAG(2)  /  32767.  -1  / 

C 

C  DATA  MCHEPS( 1 ) .MCHEPS(2)  /  0032400.  OOOOOOO  / 

C  DATA  MINMAG( 1 ) .MINMAG( 2  )  /  0000200.  OOOOOOO  / 

C  data  MAXMAG( 1 ).M4XMAG(2)  /  0077777.  0177777  / 

C  MACHINE  constants  FOR  THE  UNIVAC  1100  SERIES. 

C 

C  DATA  RMACH(  1  I  /  0147400000000  / 

C  DATA  RMACH(2)  /  0000400000000  / 

C  DATA  RMACHO)  /  0377777777777  / 

C 

C  MACHINE  constants  FOR  THE  VAX  11/780 

C  (EXPRESSED  IN  INTEGER  AND  HEXADECIMAL) 

C 

C  DATA  MCHEPS( 1 )  /  13568  / 

C  DATA  MINMAGI 1 )  /  128  / 

C  DATA  MAXMAG( 1 )  /  -32769  / 

C 

C  DATA  MCHEPS(1)  /  Z00003500  / 

C  DATA  MINMAG(1)  /  Z00000080  / 

C  DATA  MAXMAG( 1 )  /  ZFFFF7FFF  / 

C 

SPMPAR  =  RMACH(I) 

RETURN 

C 
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C  LAST  CARD  OF  FUNCTION  SPMPAR . 

C 

END 
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